Data: NIRCam simulated images obtained using MIRAGE and run through the JWST pipeline of the Large Magellanic Cloud (LMC) Astrometric Calibration Field. Simulations is obtained using a 4-pt subpixel dither for three couples of wide filters: F070W, F115W, and F200W for the SW channel, and F277W, F356W, and F444W for the LW channel. We simulated only 1 NIRCam SW detector (i.e., "NRCB1").
For this example, we use Level-2 images (.cal, calibrated but not rectified) for two SW filters (i.e., F115W and F200W) and derive the photometry in each one of them. The images for the other filters are also available and can be used to test the notebook and/or different filters combination.
The notebook is divided in two parts: in Part I we show how to create a PSF model and perform the PSF photometry, whereas in Part II, we show how to derive the final calibrated Color-Magnitude Diagram.
PSF Photometry can be obtained using:
The notebook shows:
Final plots show:
Note on pysynphot: Data files for pysynphot are distributed separately by Calibration Reference Data System. They are expected to follow a certain directory structure under the root directory, identified by the PYSYN_CDBS environment variable that must be set prior to using this package. In the example below, the root directory is arbitrarily named /my/local/dir/trds/. \ export PYSYN_CDBS=/my/local/dir/trds/ \ See documentation here for the configuration and download of the data files.
import os
import sys
import time
import numpy as np
import pandas as pd
import glob as glob
import jwst
from jwst.datamodels import ImageModel
import tarfile
import urllib.request
from astropy import wcs
from astropy import units as u
from astropy.io import fits
from astropy.visualization import (ZScaleInterval, SqrtStretch, ImageNormalize)
from astropy.visualization import simple_norm
from astropy.nddata import Cutout2D, NDData
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.table import Table, QTable
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.wcs.utils import pixel_to_skycoord
from astropy.coordinates import SkyCoord, match_coordinates_sky
from astropy.stats import sigma_clipped_stats
from photutils import CircularAperture, EPSFBuilder, find_peaks, CircularAnnulus
from photutils.detection import DAOStarFinder, IRAFStarFinder
from photutils.psf import DAOGroup, IntegratedGaussianPRF, extract_stars, IterativelySubtractedPSFPhotometry
from photutils.background import MMMBackground, MADStdBackgroundRMS
from photutils.centroids import centroid_2dg
from photutils import aperture_photometry
from ipywidgets import interact
import webbpsf
from webbpsf.utils import to_griddedpsfmodel
import pysynphot # PYSIN_CDBS must be defined in the user's environment (see note above)
%matplotlib inline
from matplotlib import style, pyplot as plt
import matplotlib.patches as patches
import matplotlib.ticker as ticker
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 30
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 30
font1 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '12'}
font2 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '20'}
We load all the images and we create a dictionary that contains all of them, divided by detectors and filters. This is useful to check which detectors and filters are available and to decide if we want to perform the photometry on all of them or only on a subset (for example, only on the SW filters).
We also create a dictionary with some useful parameters for the analysis. The dictionary contains the photometric zeropoints (from MIRAGE configuration files) and the NIRCam point spread function (PSF) FWHM, from the NIRCam Point Spread Function JDox page. The FWHM are calculated from the analysis of the expected NIRCam PSFs simulated with WebbPSF.
Note: this dictionary will be updated once the values for zeropoints and FWHM will be available for each detectors after commissioning.
Hence, we have two dictionaries:
dict_images = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}
dict_filter_short = {}
dict_filter_long = {}
ff_short = []
det_short = []
det_long = []
ff_long = []
detlist_short = []
detlist_long = []
filtlist_short = []
filtlist_long = []
if not glob.glob('./*cal*fits'):
print("Downloading images")
boxlink_images_lev2 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level2.tar.gz'
boxfile_images_lev2 = './images_level2.tar.gz'
urllib.request.urlretrieve(boxlink_images_lev2, boxfile_images_lev2)
tar = tarfile.open(boxfile_images_lev2, 'r')
tar.extractall()
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))
else:
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))
for image in images:
im = fits.open(image)
f = im[0].header['FILTER']
d = im[0].header['DETECTOR']
if d == 'NRCBLONG':
d = 'NRCB5'
elif d == 'NRCALONG':
d = 'NRCA5'
else:
d = d
wv = np.float(f[1:3])
if wv > 24:
ff_long.append(f)
det_long.append(d)
else:
ff_short.append(f)
det_short.append(d)
detlist_short = sorted(list(dict.fromkeys(det_short)))
detlist_long = sorted(list(dict.fromkeys(det_long)))
unique_list_filters_short = []
unique_list_filters_long = []
for x in ff_short:
if x not in unique_list_filters_short:
dict_filter_short.setdefault(x, {})
for x in ff_long:
if x not in unique_list_filters_long:
dict_filter_long.setdefault(x, {})
for d_s in detlist_short:
dict_images[d_s] = dict_filter_short
for d_l in detlist_long:
dict_images[d_l] = dict_filter_long
filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))
if len(dict_images[d][f]) == 0:
dict_images[d][f] = {'images': [image]}
else:
dict_images[d][f]['images'].append(image)
print("Available Detectors for SW channel:", detlist_short)
print("Available Detectors for LW channel:", detlist_long)
print("Available SW Filters:", filtlist_short)
print("Available LW Filters:", filtlist_long)
filters = ['F070W', 'F090W', 'F115W', 'F140M', 'F150W2', 'F150W', 'F162M', 'F164N', 'F182M',
'F187N', 'F200W', 'F210M', 'F212N', 'F250M', 'F277W', 'F300M', 'F322W2', 'F323N',
'F335M', 'F356W', 'F360M', 'F405N', 'F410M', 'F430M', 'F444W', 'F460M', 'F466N', 'F470N', 'F480M']
psf_fwhm = [0.987, 1.103, 1.298, 1.553, 1.628, 1.770, 1.801, 1.494, 1.990, 2.060, 2.141, 2.304, 2.341, 1.340,
1.444, 1.585, 1.547, 1.711, 1.760, 1.830, 1.901, 2.165, 2.179, 2.300, 2.302, 2.459, 2.507, 2.535, 2.574]
zp_modA = [25.7977, 25.9686, 25.8419, 24.8878, 27.0048, 25.6536, 24.6957, 22.3073, 24.8258, 22.1775, 25.3677, 24.3296,
22.1036, 22.7850, 23.5964, 24.8239, 23.6452, 25.3648, 20.8604, 23.5873, 24.3778, 23.4778, 20.5588,
23.2749, 22.3584, 23.9731, 21.9502, 20.0428, 19.8869, 21.9002]
zp_modB = [25.7568, 25.9771, 25.8041, 24.8738, 26.9821, 25.6279, 24.6767, 22.2903, 24.8042, 22.1499, 25.3391, 24.2909,
22.0574, 22.7596, 23.5011, 24.6792, 23.5769, 25.3455, 20.8631, 23.4885, 24.3883, 23.4555, 20.7007,
23.2763, 22.4677, 24.1562, 22.0422, 20.1430, 20.0173, 22.4086]
dict_utils = {filters[i]: {'psf fwhm': psf_fwhm[i], 'VegaMAG zp modA': zp_modA[i],
'VegaMAG zp modB': zp_modB[i]} for i in range(len(filters))}
If we are interested only in some filters (and/or some detectors) in the analysis, as in this example, we can select the Level-2 calibrated images from the dictionary for those filters (detectors) and analyze only those images.
In this particular example, we analyze images for filters F115W and F200W for the detector NRCB1.
dets_short = ['NRCB1'] # detector of interest in this example
filts_short = ['F115W', 'F200W'] # filters of interest in this example
To check that our images do not present artifacts and can be used in the analysis, we display them using an interactive cursor that allows to shuffle through the different images for each filter.
this is only a sketch of what I would like to show (I am not very familiar with ipywidgets). Would it be possible to show both filters at the same time, in a 2 window panel as in the static plot below? Or even better, have a widget control that allows to select the filters available and then use interact to cycle through the images?
# cell for display images using ipywidgets
def browse_images(images):
n = len(images)
def view_image(image):
det = 'NRCB1'
filt = 'F115W'
im = fits.open(dict_images[det][filt]['images'][image])
data_sb = im[1].data
norm = simple_norm(data_sb, 'sqrt', percent=99.)
plt.figure(figsize=(10, 10))
plt.title(filt)
plt.imshow(data_sb, norm=norm, cmap='Greys')
plt.show()
interact(view_image, image=(0, n - 1))
browse_images(dict_images['NRCB1']['F115W']['images'])
Cell below should be removed once we finalize the interactive one above.
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(1, len(filts_short), i + 1)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(data_sb, norm=norm, cmap='Greys')
plt.tight_layout()
We create a dictionary that contains the PSF created using WebbPSF for the detectors and filters selected above.
dict_psfs_webbpsf = {}
for det in dets_short:
dict_psfs_webbpsf.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_psfs_webbpsf[det].setdefault(filt, {})
dict_psfs_webbpsf[det][filt]['psf model grid'] = None
dict_psfs_webbpsf[det][filt]['psf model single'] = None
The function below allows to create a single PSF or a grid of PSFs and allows to save the PSF as a fits file. The model PSF are stored by default in the psf dictionary. For the grid of PSFs, users can select the number of PSFs to be created. The PSF can be created detector sampled or oversampled (the oversample can be changed inside the function).
Note: The default source spectrum is, if pysynphot is installed, a G2V star spectrum from Castelli & Kurucz (2004). Without pysynphot, the default is a simple flat spectrum such that the same number of photons are detected at each wavelength.
def create_psf_model(fov=11, create_grid=False, num=9, save_psf=False, detsampled=False):
nrc = webbpsf.NIRCam()
nrc.detector = det
nrc.filter = filt
src = webbpsf.specFromSpectralType('G5V', catalog='phoenix')
if detsampled:
print("Creating a detector sampled PSF")
aa = 'detector sampled'
fov = 21
else:
print("Creating a oversampled PSF")
aa = 'oversampled'
fov = fov
print("Using a {field}".format(field=fov), "px fov")
if create_grid:
print("")
print("Creating a grid of PSF for filter {filt} and detector {det}".format(filt=filt, det=det))
print("")
num = num
if save_psf:
outname = "./PSF_%s_samp4_G5V_fov%d_npsfs%d.fits" % (filt, fov, num)
nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,
save=True, outfile=outname, use_detsampled_psf=detsampled)
else:
grid_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,
fov_pixels=fov, use_detsampled_psf=detsampled)
dict_psfs_webbpsf[det][filt]['psf model grid'] = grid_psf
else:
print("")
print("Creating a single PSF for filter {filt} and detector {det}".format(filt=filt, det=det))
print("")
num = 1
if save_psf:
outname = "./PSF_%s_samp4_G5V_fov%d_npsfs%d.fits" % (filt, fov, num)
nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,
save=True, outfile=outname, use_detsampled_psf=detsampled)
else:
single_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,
fov_pixels=fov, use_detsampled_psf=detsampled)
dict_psfs_webbpsf[det][filt]['psf model single'] = single_psf
return
for det in dets_short:
for filt in filts_short:
create_psf_model(fov=11, num=25, create_grid=False, save_psf=False, detsampled=False)
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
ax = plt.subplot(1, 2, i + 1)
norm_epsf = simple_norm(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], 'log', percent=99.)
ax.set_title(filt, fontsize=40)
ax.imshow(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], norm=norm_epsf)
ax.set_xlabel('X [px]', fontsize=30)
ax.set_ylabel('Y [px]', fontsize=30)
plt.tight_layout()
for det in dets_short:
for filt in filts_short:
create_psf_model(fov=11, num=25, create_grid=True, save_psf=False, detsampled=False)
We show for 1 filter (F115W) the grid of PSFs and the difference from the mean
webbpsf.gridded_library.display_psf_grid(dict_psfs_webbpsf[dets_short[0]][filts_short[0]]['psf model grid'],
zoom_in=False, figsize=(14, 14))
More information on the PhotUtils Effective PSF can be found here.
threshold (approximately; threshold is applied to a convolved image) and have a size and shape similar to the defined 2D Gaussian kernel. \
Note: The threshold and the maximum distance to the closest neighbour depend on the user science case (i.e.; number of stars in the field of view, crowding, number of bright sources, minimum number of stars required to build the ePSF, etc.) and must be modified accordingly. We create a dictionary that contains the effective PSF for the detectors and filters selected above.
dict_psfs_epsf = {}
for det in dets_short:
dict_psfs_epsf.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_psfs_epsf[det].setdefault(filt, {})
dict_psfs_epsf[det][filt]['table psf stars'] = {}
dict_psfs_epsf[det][filt]['epsf single'] = {}
dict_psfs_epsf[det][filt]['epsf grid'] = {}
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = None
dict_psfs_epsf[det][filt]['epsf single'][i + 1] = None
dict_psfs_epsf[det][filt]['epsf grid'][i + 1] = None
Note that the unit of the Level-2 and Level-3 Images from the pipeline is MJy/sr (hence a surface brightness). The actual unit of the image can be checked from the header keyword BUNIT. The scalar conversion constant is copied to the header keyword PHOTMJSR, which gives the conversion from DN/s to megaJy/steradian. For our analysis we revert back to DN/s.
def find_stars_epsf(det='NRCA1', filt='F070W', dist_sel=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
print("Finding PSF stars on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px")
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th[j] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
psf_stars = daofind(data)
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = psf_stars
if dist_sel:
print("")
print("Calculating closest neigbhour distance")
d = []
daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
stars_tot = daofind_tot(data)
x_tot = stars_tot['xcentroid']
y_tot = stars_tot['ycentroid']
for xx, yy in zip(psf_stars['xcentroid'], psf_stars['ycentroid']):
sep = []
dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)
sep = np.sort(dist)[1:2][0]
d.append(sep)
psf_stars['min distance'] = d
mask_dist = (psf_stars['min distance'] > min_sep[j])
psf_stars = psf_stars[mask_dist]
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = psf_stars
print("Minimum distance required:", min_sep[j], "px")
print("")
print("Number of isolated sources found in the image used to build ePSF for {f}:".format(f=filt), len(psf_stars))
print("-----------------------------------------------------")
print("")
else:
print("")
print("Number of sources used to build ePSF for {f}:".format(f=filt), len(psf_stars))
print("--------------------------------------------")
print("")
tic = time.perf_counter()
th = [700, 500] # threshold level for the two filters (length must match number of filters analyzed)
min_sep = [10, 10] # minimum separation acceptable for ePSF stars from closest neighbour
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
find_stars_epsf(det=det, filt=filt, dist_sel=False)
toc = time.perf_counter()
print("Elapsed Time for finding stars:", toc - tic)
def build_epsf(det='NRCA1', filt='F070W'):
mmm_bkg = MMMBackground()
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
data = data_sb / imh['PHOTMJSR']
hsize = (sizes[j] - 1) / 2
x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']
y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']
mask = ((x > hsize) & (x < (data.shape[1] - 1 - hsize)) & (y > hsize) & (y < (data.shape[0] - 1 - hsize)))
stars_tbl = Table()
stars_tbl['x'] = x[mask]
stars_tbl['y'] = y[mask]
bkg = mmm_bkg(data)
data_bkgsub = data.copy()
data_bkgsub -= bkg
nddata = NDData(data=data_bkgsub)
stars = extract_stars(nddata, stars_tbl, size=sizes[j])
print("Creating ePSF for image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
epsf_builder = EPSFBuilder(oversampling=oversample, maxiters=3, progress_bar=False)
epsf, fitted_stars = epsf_builder(stars)
dict_psfs_epsf[det][filt]['epsf single'][i + 1] = epsf
Note: here we limit the maximum number of iterations to 3 (to limit it’s run time), but in practice one should use about 10 or more iterations.
tic = time.perf_counter()
sizes = [11, 11] # size of the cutout (extract region) for each PSF star - must match number of filters analyzed
oversample = 4
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
build_epsf(det=det, filt=filt)
toc = time.perf_counter()
print("Time to build the Effective PSF:", toc - tic)
We display only 1 ePSF for each filter
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
ax = plt.subplot(1, 2, i + 1)
norm_epsf = simple_norm(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, 'log', percent=99.)
plt.title(filt, fontsize=30)
ax.imshow(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, norm=norm_epsf)
Two functions:
The purpose of the first function is to count how many good PSF stars are in each sub-region defined by the grid number N. The function should start from the number provided by the user and iterate until the minimum grid size 2x2. Depending on the number of PSF stars that the users want in each cell of the grid, they can choose the appropriate grid size or modify the threshold values for the stars detection, selected when creating the single ePSF (in the Finding stars cell above).
The second function creates a grid of PSFs with EPSFBuilder. The function will return a a GriddedEPSFModel object containing a 3D array of N × n × n. The 3D array represents the N number of 2D n × n ePSFs created. It should include a grid_xypos key which will state the position of the PSF on the detector for each of the PSFs. The order of the tuples in grid_xypos refers to the number the PSF is in the 3D array.
def count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40):
num_grid_calc = np.arange(2, grid_points + 1, 1)
num_grid_calc = num_grid_calc[::-1]
for num in num_grid_calc:
print("Calculating the number of PSF stars in a %d x %d grid:" % (num, num))
print("")
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
points = np.int16((data_sb.shape[0] / num) / 2)
x_center = np.arange(points, 2 * points * (num), 2 * points)
y_center = np.arange(points, 2 * points * (num), 2 * points)
centers = np.array(np.meshgrid(x_center, y_center)).T.reshape(-1, 2)
for n, val in enumerate(centers):
x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']
y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']
flux = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['flux']
half_size = (size - 1) / 2
lim1 = val[0] - points + half_size
lim2 = val[0] + points - half_size
lim3 = val[1] - points + half_size
lim4 = val[1] + points - half_size
test = (x > lim1) & (x < lim2) & (y > lim3) & (y < lim4)
# if np.count_nonzero(test) < min_numpsf:
# raise ValueError("Not enough PSF stars in all the cells (> %d): Decrease your grid size or the minimum number of PSF stars in each cell or change parameters in the finder" %(min_numpsf))
if np.count_nonzero(test) < min_numpsf:
print("Center Coordinates of grid cell %d are (%d, %d) --- Not enough PSF stars in the cell (number of PSF stars < %d)" % (i + 1, val[0], val[1], min_numpsf))
else:
print("Center Coordinate of grid cell %d are (%d, %d) --- Number of PSF stars:" % (n + 1, val[0], val[1]), np.count_nonzero(test))
print("")
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
print("Analyzing image {number} of filter {f}, detector {d} ".format(number=i + 1, f=filt, d=det))
print("")
count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40)
Here goes the function that creates a grid of ePSF that can be saved in the epsf dictionary.
Here goes the function that create a grid of PSF models obtained perturbating the WebbPSF PSF models using the ePSF grid created above.
We perform the PSF photometry on the images, saving by default the output catalogs and the residual images in the dictionary created below. It is also possible to save the output catalogs (pickles pandas object) and residual images (fits files) in the current directory using the parameters save_output and save_residuals.
dict_phot = {}
for det in dets_short:
dict_phot.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_phot[det].setdefault(filt, {})
dict_phot[det][filt]['residual images'] = {}
dict_phot[det][filt]['output photometry tables'] = {}
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
dict_phot[det][filt]['residual images'][i + 1] = None
dict_phot[det][filt]['output photometry tables'][i + 1] = None
Note: since performing the PSF photometry on the images takes some time (for the 8 images in this example ~ 4 hours), to speed up the notebook, we use a high threshold in the finding algorithm (threshold ~ 2000) and we will use in the analyis below the catalogs obtained with a sigma threshold = 10 from a previous reduction run. To perform a meaningful data reduction, the user should modify the threshold accordingly.
Here we use as PSF model the grid of WebbPSF PSFs, but the users can change the model and use the others available (i.e., single WebbPSF PSF, single ePSF) modifying the psf parameter in the function.
def psf_phot(det='NRCA1', filt='F070W', th=2000, psf='grid_webbpsf', save_residuals=False, save_output=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
fitter = LevMarLSQFitter()
im = fits.open(dict_images[det][filt]['images'][i])
imh = im[1].header
data_sb = im[1].data
d = im[0].header['DETECTOR']
prim_dith_pos = im[0].header['PATT_NUM']
prim_dith_num = im[0].header['NUMDTHPT']
subpx_dith_pos = im[0].header['SUBPXNUM']
subpx_dith_num = im[0].header['SUBPXPNS']
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
print("Applying conversion to the data")
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf)
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
daogroup = DAOGroup(5.0 * sigma_psf)
# grid PSF
if psf == 'grid_webbpsf':
print("Using as PSF model WebbPSF PSFs grid")
psf_model = dict_psfs_webbpsf[det][filt]['psf model grid'].copy()
# single psf:
if psf == 'single_webbpsf':
print("Using as PSF model WebbPSF single PSF")
psf_model = dict_psfs_webbpsf[det][filt]['psf model single'].copy()
# epsf:
if psf == 'single_epsf':
print("Using as PSF model single ePSF")
psf_model = dict_psfs_epsf[det][filt]['epsf single'][i + 1].copy()
print("Performing the photometry on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
tic = time.perf_counter()
phot = IterativelySubtractedPSFPhotometry(finder=daofind, group_maker=daogroup,
bkg_estimator=mmm_bkg, psf_model=psf_model,
fitter=LevMarLSQFitter(),
niters=2, fitshape=(11, 11), aperture_radius=ap_radius[j])
result = phot(data)
toc = time.perf_counter()
print("Time needed to perform photometry on image {number}:".format(number=i + 1), "%.2f" % ((toc - tic) / 3600), "hours")
print("Number of sources detected in image {number} for filter {f}:".format(number=i + 1, f=filt), len(result))
residual_image = phot.get_residual_image()
dict_phot[det][filt]['residual images'][i + 1] = residual_image
dict_phot[det][filt]['output photometry tables'][i + 1] = result
# save the residual images as fits file:
if save_residuals:
hdu = fits.PrimaryHDU(residual_image)
hdul = fits.HDUList([hdu])
residual_outname = 'residual_%s_%s_webbPSF_gridPSF_%dof%d_%dof%d.fits' % (d, filt, prim_dith_pos, prim_dith_num, subpx_dith_pos, subpx_dith_num)
dir_output_phot = './'
hdul.writeto(os.path.join(dir_output_phot, residual_outname))
outname = 'phot_%s_%s_webbPSF_gridPSF_level2_%dof%d_%dof%d.pkl' % (d, filt, prim_dith_pos, prim_dith_num, subpx_dith_pos, subpx_dith_num)
# save the output photometry Tables
if save_output:
tab = result.to_pandas()
tab.to_pickle(os.path.join(dir_output_phot, outname))
tic_tot = time.perf_counter()
ap_radius = [3.0, 3.5] # must match the number of filters analyzed
if glob.glob('./*residual*.fits'):
print("Deleting Residual images from directory")
files = glob.glob('./residual*.fits')
for file in files:
os.remove(file)
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
psf_phot(det=det, filt=filt, th=2000, psf='grid_webbpsf', save_residuals=True, save_output=False)
toc_tot = time.perf_counter()
print("Time elapsed to perform the photometry of the {number} images:".format(number=(len(filts_short) * len(dict_images[det][filt]['images']))), "%.2f" % ((toc_tot - tic_tot) / 3600), "hours")
dict_phot['NRCB1']['F115W']['output photometry tables'][1]
As an example, we show the comparison between one science image and the residual image after the data reduction for both filters. Note that the residual image is obtained from the photometry run in the cell above with a very high detection threshold.
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(2, len(filts_short), i + 1)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(data_sb, norm=norm, cmap='Greys')
for det in dets_short:
for i, filt in enumerate(filts_short):
res = dict_phot[det][filt]['residual images'][1]
ax = plt.subplot(2, len(filts_short), i + 3)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(res, norm=norm, cmap='Greys')
plt.tight_layout()
Note: here we use the reduction obtained using a grid of WebbPSF PSFs as PSF models. The users can perform the data analysis using different PSF models (single PSF model, PSF grid, etc.) and compare the results.
if not glob.glob('./*phot*gridPSF*.pkl'):
print("Downloading Photometry Output")
boxlink_cat_f115w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F115W.tar.gz'
boxfile_cat_f115w = './phot_cat_F115W.tar.gz'
urllib.request.urlretrieve(boxlink_cat_f115w, boxfile_cat_f115w)
tar = tarfile.open(boxfile_cat_f115w, 'r')
tar.extractall()
boxlink_cat_f200w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F200W.tar.gz'
boxfile_cat_f200w = './phot_cat_F200W.tar.gz'
urllib.request.urlretrieve(boxlink_cat_f200w, boxfile_cat_f200w)
tar = tarfile.open(boxfile_cat_f200w, 'r')
tar.extractall()
cat_dir = './'
phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))
else:
cat_dir = './'
phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))
results_f115w = []
results_f200w = []
for phot_pkl_f115w, phot_pkl_f200w in zip(phots_pkl_f115w, phots_pkl_f200w):
ph_f115w = pd.read_pickle(phot_pkl_f115w)
ph_f200w = pd.read_pickle(phot_pkl_f200w)
result_f115w = QTable.from_pandas(ph_f115w)
result_f200w = QTable.from_pandas(ph_f200w)
results_f115w.append(result_f115w)
results_f200w.append(result_f200w)
In order to assign the WCS coordinate and hence cross-match the images, we need to transform the images to DataModel. The coordinates are assigned during the step assign_wcs step in the JWST pipeline and allow us to cross-match the different catalogs obtained for each filter.
images_f115w = []
images_f200w = []
for i in np.arange(0, len(dict_images['NRCB1']['F115W']['images']), 1):
image_f115w = ImageModel(dict_images['NRCB1']['F115W']['images'][i])
images_f115w.append(image_f115w)
for i in np.arange(0, len(dict_images['NRCB1']['F200W']['images']), 1):
image_f200w = ImageModel(dict_images['NRCB1']['F200W']['images'][i])
images_f200w.append(image_f200w)
We cross-match the catalogs to obtain the single color-magnitude diagrams.
Stars from the two filters are associated if the distance between the matches is < 0.5 px.
results_clean_f115w = []
results_clean_f200w = []
for i in np.arange(0, len(images_f115w), 1):
mask_f115w = ((results_f115w[i]['x_fit'] > 0) & (results_f115w[i]['x_fit'] < 2048) &
(results_f115w[i]['y_fit'] > 0) & (results_f115w[i]['y_fit'] < 2048) &
(results_f115w[i]['flux_fit'] > 0))
result_clean_f115w = results_f115w[i][mask_f115w]
ra_f115w, dec_f115w = images_f115w[i].meta.wcs(result_clean_f115w['x_fit'], result_clean_f115w['y_fit'])
radec_f115w = SkyCoord(ra_f115w, dec_f115w, unit='deg')
result_clean_f115w['radec'] = radec_f115w
results_clean_f115w.append(result_clean_f115w)
mask_f200w = ((results_f200w[i]['x_fit'] > 0) & (results_f200w[i]['x_fit'] < 2048) &
(results_f200w[i]['y_fit'] > 0) & (results_f200w[i]['y_fit'] < 2048) &
(results_f200w[i]['flux_fit'] > 0))
result_clean_f200w = results_f200w[i][mask_f200w]
ra_f200w, dec_f200w = images_f200w[i].meta.wcs(result_clean_f200w['x_fit'], result_clean_f200w['y_fit'])
radec_f200w = SkyCoord(ra_f200w, dec_f200w, unit='deg')
result_clean_f200w['radec'] = radec_f200w
results_clean_f200w.append(result_clean_f200w)
max_sep = 0.015 * u.arcsec
matches_phot_single = []
filt1 = 'F115W'
filt2 = 'F200W'
for res1, res2 in zip(results_clean_f115w, results_clean_f200w):
idx, d2d, _ = match_coordinates_sky(res1['radec'], res2['radec'])
sep_constraint = d2d < max_sep
match_phot_single = Table()
x_0_f115w = res1['x_0'][sep_constraint]
y_0_f115w = res1['y_0'][sep_constraint]
x_fit_f115w = res1['x_fit'][sep_constraint]
y_fit_f115w = res1['y_fit'][sep_constraint]
radec_f115w = res1['radec'][sep_constraint]
mag_f115w = (-2.5 * np.log10(res1['flux_fit']))[sep_constraint]
emag_f115w = (1.086 * (res1['flux_unc'] / res1['flux_fit']))[sep_constraint]
x_0_f200w = res2['x_0'][idx[sep_constraint]]
y_0_f200w = res2['y_0'][idx[sep_constraint]]
x_fit_f200w = res2['x_fit'][idx[sep_constraint]]
y_fit_f200w = res2['y_fit'][idx[sep_constraint]]
radec_f200w = res2['radec'][idx][sep_constraint]
mag_f200w = (-2.5 * np.log10(res2['flux_fit']))[idx[sep_constraint]]
emag_f200w = (1.086 * (res2['flux_unc'] / res2['flux_fit']))[idx[sep_constraint]]
match_phot_single['x_0_' + filt1] = x_0_f115w
match_phot_single['y_0_' + filt1] = y_0_f115w
match_phot_single['x_fit_' + filt1] = x_fit_f115w
match_phot_single['y_fit_' + filt1] = y_fit_f115w
match_phot_single['radec_' + filt1] = radec_f115w
match_phot_single['mag_' + filt1] = mag_f115w
match_phot_single['emag_' + filt1] = emag_f115w
match_phot_single['x_0_' + filt2] = x_0_f200w
match_phot_single['y_0_' + filt2] = y_0_f200w
match_phot_single['x_fit_' + filt2] = x_fit_f200w
match_phot_single['y_fit_' + filt2] = y_fit_f200w
match_phot_single['radec_' + filt2] = radec_f200w
match_phot_single['mag_' + filt2] = mag_f200w
match_phot_single['emag_' + filt2] = emag_f200w
matches_phot_single.append(match_phot_single)
plt.figure(figsize=(12, 16))
plt.clf()
for i in np.arange(0, len(matches_phot_single), 1):
ax = plt.subplot(2, 2, i + 1)
j = str(i + 1)
xlim0 = -0.5
xlim1 = 0.8
ylim0 = -1
ylim1 = -9
ax.set_xlim(xlim0, xlim1)
ax.set_ylim(ylim0, ylim1)
ax.xaxis.set_major_locator(ticker.AutoLocator())
ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax.yaxis.set_major_locator(ticker.AutoLocator())
ax.yaxis.set_minor_locator(ticker.AutoMinorLocator())
f115w_single = matches_phot_single[i]['mag_' + filt1]
f200w_single = matches_phot_single[i]['mag_' + filt2]
ax.scatter(f115w_single - f200w_single, f115w_single, s=1, color='k')
ax.set_xlabel(filt1 + '-' + filt2, fontdict=font2)
ax.set_ylabel(filt1, fontdict=font2)
ax.text(xlim0 + 0.1, -8.65, "Image %s" % j, fontdict=font2)
plt.tight_layout()
We show the difference in the stars position derived from daofind and the psf fitting algorithm. We also show the difference $\Delta$X and $\Delta$Y as a function of the instrumental magnitudes.
plt.figure(figsize=(12, 6))
ax1 = plt.subplot(1, 2, 1)
xlim0 = -1
xlim1 = 1
ylim0 = -1
ylim1 = 1
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
x_find_f115w = results_clean_f115w[0]['x_0']
y_find_f115w = results_clean_f115w[0]['y_0']
x_psf_f115w = results_clean_f115w[0]['x_fit']
y_psf_f115w = results_clean_f115w[0]['y_fit']
delta_x_f115w = x_find_f115w - x_psf_f115w
delta_y_f115w = y_find_f115w - y_psf_f115w
_, d_x_f115w, sigma_d_x_f115w = sigma_clipped_stats(delta_x_f115w)
_, d_y_f115w, sigma_d_y_f115w = sigma_clipped_stats(delta_y_f115w)
ax1.scatter(delta_x_f115w, delta_y_f115w, s=1, color='gray')
ax1.set_xlabel('$\Delta$ X (px)', fontdict=font2)
ax1.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax1.set_title(filt1, fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 0.15, ' $\Delta$ X = %5.3f $\pm$ %5.3f' % (d_x_f115w, sigma_d_x_f115w),
color='k', fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 0.30, ' $\Delta$ Y = %5.3f $\pm$ %5.3f' % (d_y_f115w, sigma_d_y_f115w),
color='k', fontdict=font2)
ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2 = plt.subplot(1, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
x_find_f200w = results_clean_f200w[0]['x_0']
y_find_f200w = results_clean_f200w[0]['y_0']
x_psf_f200w = results_clean_f200w[0]['x_fit']
y_psf_f200w = results_clean_f200w[0]['y_fit']
delta_x_f200w = x_find_f200w - x_psf_f200w
delta_y_f200w = y_find_f200w - y_psf_f200w
_, d_x_f200w, sigma_d_x_f200w = sigma_clipped_stats(delta_x_f200w)
_, d_y_f200w, sigma_d_y_f200w = sigma_clipped_stats(delta_y_f200w)
ax2.scatter(delta_x_f200w, delta_y_f200w, s=1, color='gray')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, ' $\Delta$ X = %5.3f $\pm$ %5.3f' % (d_x_f200w, sigma_d_x_f200w),
color='k', fontdict=font2)
ax2.text(xlim0 + 0.05, ylim1 - 0.30, ' $\Delta$ Y = %5.3f $\pm$ %5.3f' % (d_y_f200w, sigma_d_y_f200w),
color='k', fontdict=font2)
ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.set_xlabel('$\Delta$ X (px)', fontdict=font2)
ax2.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax2.set_title(filt2, fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(12, 8))
ax1 = plt.subplot(2, 2, 1)
xlim0 = -9
xlim1 = -1
ylim0 = -1
ylim1 = 1
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
mag_inst_f115w = -2.5 * np.log10(results_clean_f115w[0]['flux_fit'])
ax1.scatter(mag_inst_f115w, delta_x_f115w, s=1, color='gray')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax1.set_xlabel(filt1 + '_inst', fontdict=font2)
ax1.set_ylabel('$\Delta$ X (px)', fontdict=font2)
ax2 = plt.subplot(2, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(mag_inst_f115w, delta_y_f115w, s=1, color='gray')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.set_xlabel(filt1 + '_inst', fontdict=font2)
ax2.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax3 = plt.subplot(2, 2, 3)
ax3.set_xlim(xlim0, xlim1)
ax3.set_ylim(ylim0, ylim1)
ax3.xaxis.set_major_locator(ticker.AutoLocator())
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.yaxis.set_major_locator(ticker.AutoLocator())
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())
mag_inst_f200w = -2.5 * np.log10(results_clean_f200w[0]['flux_fit'])
ax3.scatter(mag_inst_f200w, delta_x_f200w, s=1, color='gray')
ax3.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax3.set_xlabel(filt2 + '_inst', fontdict=font2)
ax3.set_ylabel('$\Delta$ X (px)', fontdict=font2)
ax4 = plt.subplot(2, 2, 4)
ax4.set_xlim(xlim0, xlim1)
ax4.set_ylim(ylim0, ylim1)
ax4.xaxis.set_major_locator(ticker.AutoLocator())
ax4.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax4.yaxis.set_major_locator(ticker.AutoLocator())
ax4.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax4.scatter(mag_inst_f200w, delta_y_f200w, s=1, color='gray')
ax4.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax4.set_xlabel(filt2 + '_inst', fontdict=font2)
ax4.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
plt.tight_layout()
To obtain a final color-magnitude diagram, we need to cross-match all the catalogs for each filters and then cross-match the derived final catalogs.
Note: this is the most conservative approach since we impose that a star must be found in all 4 catalogs.
I couldn't find an easier way to write this function, where you need to match the first two catalogs, derive a sub-catalogs with only the matches and then iterate for all the other catalogs available. We should also think on how to create a function that allows to keep the stars also if they are available in X out of Y catalogs (i.e., if for some reasons, a measure is not available in 1 of the images, but the star is well measured in the other 3, it doesn't make sense to discard the object).
def crossmatch_filter(table=None):
num = 0
num_cat = np.char.mod('%d', np.arange(1, len(table) + 1, 1))
idx_12, d2d_12, _ = match_coordinates_sky(table[num]['radec'], table[num + 1]['radec'])
sep_constraint_12 = d2d_12 < max_sep
matches_12 = Table()
matches_12['radec_' + num_cat[num]] = table[num]['radec'][sep_constraint_12]
matches_12['mag_' + num_cat[num]] = (-2.5 * np.log10(table[num]['flux_fit']))[sep_constraint_12]
matches_12['emag_' + num_cat[num]] = (1.086 * (table[num]['flux_unc'] /
table[num]['flux_fit']))[sep_constraint_12]
matches_12['radec_' + num_cat[num + 1]] = table[num + 1]['radec'][idx_12[sep_constraint_12]]
matches_12['mag_' + num_cat[num + 1]] = (-2.5 * np.log10(table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]
matches_12['emag_' + num_cat[num + 1]] = (1.086 * (table[num + 1]['flux_unc'] /
table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]
idx_123, d2d_123, _ = match_coordinates_sky(matches_12['radec_' + num_cat[num]], table[num + 2]['radec'])
sep_constraint_123 = d2d_123 < max_sep
matches_123 = Table()
matches_123['radec_' + num_cat[num]] = matches_12['radec_' + num_cat[num]][sep_constraint_123]
matches_123['mag_' + num_cat[num]] = matches_12['mag_' + num_cat[num]][sep_constraint_123]
matches_123['emag_' + num_cat[num]] = matches_12['emag_' + num_cat[num]][sep_constraint_123]
matches_123['radec_' + num_cat[num + 1]] = matches_12['radec_' + num_cat[num + 1]][sep_constraint_123]
matches_123['mag_' + num_cat[num + 1]] = matches_12['mag_' + num_cat[num + 1]][sep_constraint_123]
matches_123['emag_' + num_cat[num + 1]] = matches_12['emag_' + num_cat[num + 1]][sep_constraint_123]
matches_123['radec_' + num_cat[num + 2]] = table[num + 2]['radec'][idx_123[sep_constraint_123]]
matches_123['mag_' + num_cat[num + 2]] = (-2.5 * np.log10(table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]
matches_123['emag_' + num_cat[num + 2]] = (1.086 * (table[num + 2]['flux_unc'] /
table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]
idx_1234, d2d_1234, _ = match_coordinates_sky(matches_123['radec_' + num_cat[num]], table[num + 3]['radec'])
sep_constraint_1234 = d2d_1234 < max_sep
matches_1234 = Table()
matches_1234['radec_' + num_cat[num]] = matches_123['radec_' + num_cat[num]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num]] = matches_123['mag_' + num_cat[num]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num]] = matches_123['emag_' + num_cat[num]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 1]] = matches_123['radec_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num + 1]] = matches_123['mag_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num + 1]] = matches_123['emag_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 2]] = matches_123['radec_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num + 2]] = matches_123['mag_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num + 2]] = matches_123['emag_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 3]] = table[num + 3]['radec'][idx_1234[sep_constraint_1234]]
matches_1234['mag_' + num_cat[num + 3]] = (-2.5 * np.log10(table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]
matches_1234['emag_' + num_cat[num + 3]] = (1.086 * (table[num + 3]['flux_unc'] /
table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]
matches_1234
return matches_1234
matches_f115w = crossmatch_filter(table=results_clean_f115w)
matches_f200w = crossmatch_filter(table=results_clean_f200w)
For the final catalog, we assume that the magnitude is the mean of the 4 measures and the error on the magnitude is its standard deviation.
To easily perform this arithmetic operation on the table, we convert the table to pandas dataframe.
df_f115w = matches_f115w.to_pandas()
df_f200w = matches_f200w.to_pandas()
df_f115w['RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)
df_f115w['e_RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)
df_f115w['Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)
df_f115w['e_Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)
df_f115w[filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)
df_f115w['e' + filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)
df_f200w['RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)
df_f200w['e_RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)
df_f200w['Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)
df_f200w['e_Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)
df_f200w[filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)
df_f200w['e' + filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)
plt.figure(figsize=(12, 14))
plt.clf()
ax1 = plt.subplot(1, 2, 1)
ax1.set_xlabel(filt1 + '_inst -' + filt2 + '_inst', fontdict=font2)
ax1.set_ylabel(filt1 + '_inst', fontdict=font2)
xlim0 = -0.5
xlim1 = 0.8
ylim0 = -1.5
ylim1 = -9
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
radec_f115w_inst = SkyCoord(df_f115w['RA_' + filt1], df_f115w['Dec_' + filt1], unit='deg')
radec_f200w_inst = SkyCoord(df_f200w['RA_' + filt2], df_f200w['Dec_' + filt2], unit='deg')
idx_inst, d2d_inst, _ = match_coordinates_sky(radec_f115w_inst, radec_f200w_inst)
sep_constraint_inst = d2d_inst < max_sep
f115w_inst = np.array(df_f115w[filt1 + '_inst'][sep_constraint_inst])
ef115w_inst = np.array(df_f115w['e' + filt1 + '_inst'][sep_constraint_inst])
radec_f115w = radec_f115w_inst[sep_constraint_inst]
f200w_inst = np.array(df_f200w[filt2 + '_inst'][idx_inst[sep_constraint_inst]])
ef200w_inst = np.array(df_f200w['e' + filt2 + '_inst'][idx_inst[sep_constraint_inst]])
radec_f200w = radec_f200w_inst[idx_inst[sep_constraint_inst]]
ax1.scatter(f115w_inst - f200w_inst, f115w_inst, s=1, color='k')
ax2 = plt.subplot(2, 2, 2)
ax2.set_xlabel(filt1 + '_inst', fontdict=font2)
ax2.set_ylabel('$\sigma$' + filt1, fontdict=font2)
xlim0 = -9
xlim1 = -1.5
ylim0 = -0.01
ylim1 = 1
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(df_f115w[filt1 + '_inst'], df_f115w['e' + filt1 + '_inst'], s=1, color='k')
ax3 = plt.subplot(2, 2, 4)
ax3.set_xlabel(filt2 + '_inst', fontdict=font2)
ax3.set_ylabel('$\sigma$' + filt2, fontdict=font2)
ax3.set_xlim(xlim0, xlim1)
ax3.set_ylim(ylim0, ylim1)
ax3.xaxis.set_major_locator(ticker.AutoLocator())
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.yaxis.set_major_locator(ticker.AutoLocator())
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.scatter(df_f200w[filt2 + '_inst'], df_f200w['e' + filt2 + '_inst'], s=1, color='k')
plt.tight_layout()
To obtain the final calibrated color-magnitude diagram, we need to calculate the photometric zeropoints. Hence we need to perform aperture photometry on the calibrated images (Level-3), apply the appropriate aperture correction for the finite aperture adopted (the values provided in the dictionary above are for an infinite aperture) and then compare it with the PSF photometry. Hence, we can summarize the steps as follows:
dict_images_combined = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}
dict_filter_short = {}
dict_filter_long = {}
ff_short = []
det_short = []
det_long = []
ff_long = []
detlist_short = []
detlist_long = []
filtlist_short = []
filtlist_long = []
if not glob.glob('./*combined*fits'):
print("Downloading images")
boxlink_images_lev3 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level3.tar.gz'
boxfile_images_lev3 = './images_level3.tar.gz'
urllib.request.urlretrieve(boxlink_images_lev3, boxfile_images_lev3)
tar = tarfile.open(boxfile_images_lev3, 'r')
tar.extractall()
images_dir = './'
files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
else:
images_dir = './'
files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
for file in files_singles:
im = fits.open(file)
f = im[0].header['FILTER']
d = im[0].header['DETECTOR']
if d == 'NRCBLONG':
d = 'NRCB5'
elif d == 'NRCALONG':
d = 'NRCA5'
else:
d = d
wv = np.float(f[1:3])
if wv > 24:
ff_long.append(f)
det_long.append(d)
else:
ff_short.append(f)
det_short.append(d)
detlist_short = sorted(list(dict.fromkeys(det_short)))
detlist_long = sorted(list(dict.fromkeys(det_long)))
unique_list_filters_short = []
unique_list_filters_long = []
for x in ff_short:
if x not in unique_list_filters_short:
dict_filter_short.setdefault(x, {})
for x in ff_long:
if x not in unique_list_filters_long:
dict_filter_long.setdefault(x, {})
for d_s in detlist_short:
dict_images_combined[d_s] = dict_filter_short
for d_l in detlist_long:
dict_images_combined[d_l] = dict_filter_long
filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))
if len(dict_images_combined[d][f]) == 0:
dict_images_combined[d][f] = {'images': [file]}
else:
dict_images_combined[d][f]['images'].append(file)
print("Available Detectors for SW channel:", detlist_short)
print("Available Detectors for LW channel:", detlist_long)
print("Available SW Filters:", filtlist_short)
print("Available LW Filters:", filtlist_long)
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images_combined[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(1, len(filts_short), i + 1)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
ax.imshow(data_sb, norm=norm, cmap='Greys')
plt.tight_layout()
As we have done previously, we create a dictionary that contains the tables with the derived aperture photometry for each image.
dict_aper = {}
for det in dets_short:
dict_aper.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_aper[det].setdefault(filt, {})
dict_aper[det][filt]['stars for ap phot'] = None
dict_aper[det][filt]['stars for ap phot matched'] = None
dict_aper[det][filt]['aperture phot table'] = None
def find_bright_stars(det='NRCA1', filt='F070W', dist_sel=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
image = fits.open(dict_images_combined[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
print("Selecting stars for aperture photometry on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px")
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th[j] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
apcorr_stars = daofind(data)
dict_aper[det][filt]['stars for ap phot'] = apcorr_stars
if dist_sel:
print("")
print("Calculating closest neigbhour distance")
d = []
daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
stars_tot = daofind_tot(data)
x_tot = stars_tot['xcentroid']
y_tot = stars_tot['ycentroid']
for xx, yy in zip(apcorr_stars['xcentroid'], apcorr_stars['ycentroid']):
sep = []
dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)
sep = np.sort(dist)[1:2][0]
d.append(sep)
apcorr_stars['min distance'] = d
mask_dist = (apcorr_stars['min distance'] > min_sep[j])
apcorr_stars = apcorr_stars[mask_dist]
dict_aper[det][filt]['stars for ap phot'] = apcorr_stars
print("Minimum distance required:", min_sep[j], "px")
print("")
print("Number of bright isolated sources found in the image for {f}:".format(f=filt), len(apcorr_stars))
print("-----------------------------------------------------")
print("")
else:
print("")
print("Number of bright sources found in the image for {f}:".format(f=filt), len(apcorr_stars))
print("--------------------------------------------")
print("")
return
tic = time.perf_counter()
th = [700, 500] # threshold level for the two filters (length must match number of filters analyzed)
min_sep = [10, 10] # minimum separation acceptable for zp stars from closest neighbour
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):
find_bright_stars(det=det, filt=filt, dist_sel=False)
toc = time.perf_counter()
print("Elapsed Time for finding stars for Aperture Photometry:", toc - tic)
As a further way to obtain a good quality sample, we cross-match the catalogs from the two filters and retain only the stars in common
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):
image = ImageModel(dict_images_combined[det][filt]['images'][i])
ra, dec = image.meta.wcs(dict_aper[det][filt]['stars for ap phot']['xcentroid'],
dict_aper[det][filt]['stars for ap phot']['ycentroid'])
radec = SkyCoord(ra, dec, unit='deg')
dict_aper[det][filt]['stars for ap phot']['radec'] = radec
idx_ap, d2d_ap, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot']['radec'],
dict_aper[det][filt2]['stars for ap phot']['radec'])
sep_constraint_ap = d2d_ap < max_sep
matched_apcorr_f115w = Table()
matched_apcorr_f200w = Table()
matched_apcorr_f115w = dict_aper[det][filt1]['stars for ap phot'][sep_constraint_ap]
matched_apcorr_f200w = dict_aper[det][filt2]['stars for ap phot'][idx_ap[sep_constraint_ap]]
dict_aper[det][filt1]['stars for ap phot matched'] = matched_apcorr_f115w
dict_aper[det][filt2]['stars for ap phot matched'] = matched_apcorr_f200w
Note: these values are obtained from the study of the synthetic WebbPSF PSFs. They will be updated once we have in-flight measures.
if os.path.isfile('./aperture_correction_table.txt'):
ap_tab = './aperture_correction_table.txt'
else:
print("Downloading the aperture correction table")
boxlink_apcorr_table = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/aperture_correction_table.txt'
boxfile_apcorr_table = './aperture_correction_table.txt'
urllib.request.urlretrieve(boxlink_apcorr_table, boxfile_apcorr_table)
ap_tab = './aperture_correction_table.txt'
aper_table = pd.read_csv(ap_tab, header=None, sep='\s+', index_col=0,
names=['filter', 'pupil', 'wave', 'r10', 'r20', 'r30', 'r40', 'r50', 'r60', 'r70', 'r80',
'r85', 'r90', 'sky_flux_px', 'apcorr10', 'apcorr20', 'apcorr30', 'apcorr40',
'apcorr50', 'apcorr60', 'apcorr70', 'apcorr80', 'apcorr85', 'apcorr90', 'sky_in',
'sky_out'], comment='#', skiprows=0, usecols=range(0, 26))
aper_table.head()
def aperture_phot(det=det, filt='F070W'):
radii = [aper_table.loc[filt]['r70']]
ees = '70'.split()
ee_radii = dict(zip(ees, radii))
positions = np.transpose((dict_aper[det][filt]['stars for ap phot matched']['xcentroid'],
dict_aper[det][filt]['stars for ap phot matched']['ycentroid']))
image = fits.open(dict_images_combined[det][filt]['images'][0])
data_sb = image[1].data
imh = image[1].header
data = data_sb / imh['PHOTMJSR']
# sky from the aperture correction table:
sky = {"sky_in": aper_table.loc[filt]['r80'], "sky_out": aper_table.loc[filt]['r85']}
tic = time.perf_counter()
table_aper = Table()
for ee, radius in ee_radii.items():
print("Performing aperture photometry for radius equivalent to EE = {0}% for filter {1}".format(ee, filt))
aperture = CircularAperture(positions, r=radius)
annulus_aperture = CircularAnnulus(positions, r_in=sky["sky_in"], r_out=sky["sky_out"])
annulus_mask = annulus_aperture.to_mask(method='center')
bkg_median = []
for mask in annulus_mask:
annulus_data = mask.multiply(data)
annulus_data_1d = annulus_data[mask.data > 0]
_, median_sigclip, _ = sigma_clipped_stats(annulus_data_1d)
bkg_median.append(median_sigclip)
bkg_median = np.array(bkg_median)
phot = aperture_photometry(data, aperture, method='exact')
phot['annulus_median'] = bkg_median
phot['aper_bkg'] = bkg_median * aperture.area
phot['aper_sum_bkgsub'] = phot['aperture_sum'] - phot['aper_bkg']
apcorr = [aper_table.loc[filt]['apcorr70']]
phot['aper_sum_corrected'] = phot['aper_sum_bkgsub'] * apcorr
phot['mag_corrected'] = -2.5 * np.log10(phot['aper_sum_corrected']) + dict_utils[filt]['VegaMAG zp modB']
table_aper.add_column(phot['aperture_sum'], name='aper_sum_' + ee)
table_aper.add_column(phot['annulus_median'], name='annulus_median_' + ee)
table_aper.add_column(phot['aper_bkg'], name='aper_bkg_ee_' + ee)
table_aper.add_column(phot['aper_sum_bkgsub'], name='aper_sum_bkgsub_' + ee)
table_aper.add_column(phot['aper_sum_corrected'], name='aper_sum_corrected_' + filt)
table_aper.add_column(phot['mag_corrected'], name='mag_corrected_' + filt)
dict_aper[det][filt]['aperture phot table'] = table_aper
toc = time.perf_counter()
print("Time Elapsed:", toc - tic)
return
aperture_phot(det=det, filt=filt1)
aperture_phot(det=det, filt=filt2)
plt.figure(figsize=(14, 8))
plt.clf()
ax1 = plt.subplot(2, 1, 1)
ax1.set_xlabel(filt1, fontdict=font2)
ax1.set_ylabel('Zeropoint', fontdict=font2)
idx_zp_1, d2d_zp_1, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot matched']['radec'], radec_f115w_inst)
sep_constraint_zp_1 = d2d_zp_1 < max_sep
f115w_ap_matched = np.array(dict_aper[det][filt1]['aperture phot table']['mag_corrected_' + filt1][sep_constraint_zp_1])
f115w_psf_matched = np.array(df_f115w[filt1 + '_inst'][idx_zp_1[sep_constraint_zp_1]])
diff_f115w = f115w_ap_matched - f115w_psf_matched
_, zp_f115w, zp_sigma_f115w = sigma_clipped_stats(diff_f115w)
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f115w) - 0.5
ylim1 = np.mean(diff_f115w) + 0.5
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(f115w_psf_matched, diff_f115w, s=50, color='k')
ax1.plot([xlim0, xlim1], [zp_f115w, zp_f115w], color='r', lw=5, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + ' Zeropoint = %5.3f $\pm$ %5.3f' % (zp_f115w, zp_sigma_f115w), color='k', fontdict=font2)
ax2 = plt.subplot(2, 1, 2)
ax2.set_xlabel(filt2, fontdict=font2)
ax2.set_ylabel('Zeropoint', fontdict=font2)
idx_zp_2, d2d_zp_2, _ = match_coordinates_sky(dict_aper[det][filt2]['stars for ap phot matched']['radec'], radec_f200w_inst)
sep_constraint_zp_2 = d2d_zp_2 < max_sep
f200w_ap_matched = np.array(dict_aper[det][filt2]['aperture phot table']['mag_corrected_' + filt2][sep_constraint_zp_2])
f200w_psf_matched = np.array(df_f200w[filt2 + '_inst'][idx_zp_2[sep_constraint_zp_2]])
diff_f200w = f200w_ap_matched - f200w_psf_matched
_, zp_f200w, zp_sigma_f200w = sigma_clipped_stats(diff_f200w)
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f200w) - 0.5
ylim1 = np.mean(diff_f200w) + 0.5
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(f200w_psf_matched, diff_f200w, s=50, color='k')
ax2.plot([xlim0, xlim1], [zp_f200w, zp_f200w], color='r', lw=5, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + ' Zeropoint = %5.3f $\pm$ %5.3f' % (zp_f200w, zp_sigma_f200w), color='k', fontdict=font2)
plt.tight_layout()
if os.path.isfile('./pointsource.cat'):
input_cat = './pointsource.cat'
else:
print("Downloading input pointsource catalog")
boxlink_input_cat = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/pointsource.cat'
boxfile_input_cat = './pointsource.cat'
urllib.request.urlretrieve(boxlink_input_cat, boxfile_input_cat)
input_cat = './pointsource.cat'
cat = pd.read_csv(input_cat, header=None, sep='\s+', names=['ra_in', 'dec_in', 'f070w_in', 'f115w_in',
'f200w_in', 'f277w_in', 'f356w_in', 'f444w_in'],
comment='#', skiprows=7, usecols=range(0, 8))
cat.head()
Extract from the input catalog the stars in the same region as the one analyzed
lim_ra_min = np.min(radec_f115w.ra)
lim_ra_max = np.max(radec_f115w.ra)
lim_dec_min = np.min(radec_f115w.dec)
lim_dec_max = np.max(radec_f115w.dec)
cat_sel = cat[(cat['ra_in'] > lim_ra_min) & (cat['ra_in'] < lim_ra_max) & (cat['dec_in'] > lim_dec_min)
& (cat['dec_in'] < lim_dec_max)]
plt.figure(figsize=(12, 14))
plt.clf()
ax1 = plt.subplot(1, 2, 1)
mag1_in = np.array(cat_sel['f115w_in'])
mag2_in = np.array(cat_sel['f200w_in'])
diff_in = mag1_in - mag2_in
xlim0 = -0.25
xlim1 = 1.75
ylim0 = 25
ylim1 = 15
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(mag1_in - mag2_in, mag1_in, s=1, color='k')
ax1.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)
ax1.set_ylabel(filt1, fontdict=font2)
ax1.text(xlim0 + 0.15, 15.5, "Input", fontdict=font2)
ax2 = plt.subplot(1, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
f115w = f115w_inst + zp_f115w
f200w = f200w_inst + zp_f200w
maglim = np.arange(18, 25, 1)
mags = []
errs_mag = []
errs_col = []
for i in np.arange(0, len(maglim) - 1, 1):
mag = (maglim[i] + maglim[i + 1]) / 2
err_mag1 = ef115w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]
err_mag2 = ef200w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]
err_mag = np.mean(err_mag1[i])
err_temp = np.sqrt(err_mag1**2 + err_mag2**2)
err_col = np.mean(err_temp[i])
errs_mag.append(err_mag)
errs_col.append(err_col)
mags.append(mag)
col = [0] * (len(maglim) - 1)
ax2.errorbar(col, mags, yerr=errs_mag, xerr=errs_col, fmt='o', color='k')
ax2.scatter(f115w - f200w, f115w, s=1, color='k')
ax2.text(xlim0 + 0.15, 15.5, "Output", fontdict=font2)
ax2.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)
ax2.set_ylabel(filt1, fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(14, 8))
plt.clf()
ax1 = plt.subplot(2, 1, 1)
ax1.set_xlabel(filt1, fontdict=font2)
ax1.set_ylabel('$\Delta$ Mag', fontdict=font2)
radec_input = SkyCoord(cat_sel['ra_in'], cat_sel['dec_in'], unit='deg')
idx_f115w_cfr, d2d_f115w_cfr, _ = match_coordinates_sky(radec_input, radec_f115w)
sep_f115w_cfr = d2d_f115w_cfr < max_sep
f115w_inp_cfr = np.array(cat_sel['f115w_in'][sep_f115w_cfr])
f115w_psf_cfr = np.array(f115w[idx_f115w_cfr[sep_f115w_cfr]])
diff_f115w_cfr = f115w_inp_cfr - f115w_psf_cfr
_, med_diff_f115w_cfr, sig_diff_f115w_cfr = sigma_clipped_stats(diff_f115w_cfr)
xlim0 = 16
xlim1 = 24.5
ylim0 = np.mean(diff_f115w_cfr) - 0.5
ylim1 = np.mean(diff_f115w_cfr) + 0.5
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(f115w_psf_cfr, diff_f115w_cfr, s=5, color='k')
ax1.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f'
% (med_diff_f115w_cfr, sig_diff_f115w_cfr), color='k', fontdict=font2)
ax2 = plt.subplot(2, 1, 2)
ax2.set_xlabel(filt2, fontdict=font2)
ax2.set_ylabel('$\Delta$ Mag', fontdict=font2)
idx_f200w_cfr, d2d_f200w_cfr, _ = match_coordinates_sky(radec_input, radec_f200w)
sep_f200w_cfr = d2d_f200w_cfr < max_sep
f200w_inp_cfr = np.array(cat_sel['f200w_in'][sep_f200w_cfr])
f200w_psf_cfr = np.array(f200w[idx_f200w_cfr[sep_f200w_cfr]])
diff_f200w_cfr = f200w_inp_cfr - f200w_psf_cfr
_, med_diff_f200w_cfr, sig_diff_f200w_cfr = sigma_clipped_stats(diff_f200w_cfr)
xlim0 = 16
xlim1 = 24
ylim0 = np.mean(diff_f200w_cfr) - 0.5
ylim1 = np.mean(diff_f200w_cfr) + 0.5
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(f200w_psf_cfr, diff_f200w_cfr, s=5, color='k')
ax2.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f'
% (med_diff_f200w_cfr, sig_diff_f200w_cfr), color='k', fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(12, 6))
ax1 = plt.subplot(1, 2, 1)
xlim0 = -10
xlim1 = 10
ylim0 = -10
ylim1 = 10
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.set_xlabel('$\Delta$ RA (mas)', fontdict=font2)
ax1.set_ylabel('$\Delta$ Dec (mas)', fontdict=font2)
ax1.set_title(filt1, fontdict=font2)
ra_f115w_inp_cfr = np.array(cat_sel['ra_in'][sep_f115w_cfr])
ra_f115w_psf_cfr = np.array(radec_f115w.ra[idx_f115w_cfr[sep_f115w_cfr]])
dec_f115w_inp_cfr = np.array(cat_sel['dec_in'][sep_f115w_cfr])
dec_f115w_psf_cfr = np.array(radec_f115w.dec[idx_f115w_cfr[sep_f115w_cfr]])
dec_rad_f115w = np.radians(dec_f115w_psf_cfr)
diffra_f115w_cfr = ((((ra_f115w_inp_cfr - ra_f115w_psf_cfr) * np.cos(dec_rad_f115w)) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffra_f115w_cfr, sig_diffra_f115w_cfr = sigma_clipped_stats(diffra_f115w_cfr)
diffdec_f115w_cfr = (((dec_f115w_inp_cfr - dec_f115w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffdec_f115w_cfr, sig_diffdec_f115w_cfr = sigma_clipped_stats(diffdec_f115w_cfr)
ax1.scatter(diffra_f115w_cfr, diffdec_f115w_cfr, s=1, color='k')
ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 1.50, ' $\Delta$ RA (mas) = %5.3f $\pm$ %5.3f'
% (med_diffra_f115w_cfr, sig_diffra_f115w_cfr), color='k', fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 3.0, ' $\Delta$ Dec (mas) = %5.3f $\pm$ %5.3f'
% (med_diffdec_f115w_cfr, sig_diffdec_f115w_cfr), color='k', fontdict=font2)
ax2 = plt.subplot(1, 2, 2)
xlim0 = -10
xlim1 = 10
ylim0 = -10
ylim1 = 10
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.set_title(filt2, fontdict=font2)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.set_xlabel('$\Delta$ RA (mas)', fontdict=font2)
ax2.set_ylabel('$\Delta$ Dec (mas)', fontdict=font2)
ra_f200w_inp_cfr = np.array(cat_sel['ra_in'][sep_f200w_cfr])
ra_f200w_psf_cfr = np.array(radec_f200w.ra[idx_f200w_cfr[sep_f200w_cfr]])
dec_f200w_inp_cfr = np.array(cat_sel['dec_in'][sep_f200w_cfr])
dec_f200w_psf_cfr = np.array(radec_f200w.dec[idx_f200w_cfr[sep_f200w_cfr]])
dec_rad_f200w = np.radians(dec_f200w_psf_cfr)
diffra_f200w_cfr = ((((ra_f200w_inp_cfr - ra_f200w_psf_cfr) * np.cos(dec_rad_f200w)) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffra_f200w_cfr, sig_diffra_f200w_cfr = sigma_clipped_stats(diffra_f200w_cfr)
diffdec_f200w_cfr = (((dec_f200w_inp_cfr - dec_f200w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffdec_f200w_cfr, sig_diffdec_f200w_cfr = sigma_clipped_stats(diffdec_f200w_cfr)
ax2.scatter(diffra_f200w_cfr, diffdec_f200w_cfr, s=1, color='k')
ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 1.50, ' $\Delta$ RA (mas) = %5.3f $\pm$ %5.3f'
% (med_diffra_f200w_cfr, sig_diffra_f200w_cfr), color='k', fontdict=font2)
ax2.text(xlim0 + 0.05, ylim1 - 3.0, ' $\Delta$ Dec (mas) = %5.3f $\pm$ %5.3f'
% (med_diffdec_f200w_cfr, sig_diffdec_f200w_cfr), color='k', fontdict=font2)
plt.tight_layout()
This notebook provides a general overview on how to perform PSF photometry using the PhotUtils package. The choice of the different parameters adopted in all the reduction steps as well as the choice of the PSF model depend on the specific user science case. Moreover, a detailed analysis that allow to provide recommendations on how to set those parameters and outline the differences in the output photometry when different PSF models are adopted (single vs PSF grid, number of PSFs in the grid, etc.) will be possible only when real data will be available after the instrument commissioning. In this context, we note that one of the selected ERS program (ERS 1334 - The Resolved Stellar Populations Early Release Science Program) will provide a fundamental test benchmark to explore how the different choices outlined above will impact the quality of the PSF photometry in a crowded stellar region.
Author: Matteo Correnti, JWST/NIRCam STScI Scientist II \ Updated on: 2021-01-15
import os
import sys
import time
import numpy as np
import pandas as pd
import glob as glob
import jwst
from jwst.datamodels import ImageModel
import tarfile
import urllib.request
from astropy import wcs
from astropy import units as u
from astropy.io import fits
from astropy.visualization import (ZScaleInterval, SqrtStretch, ImageNormalize)
from astropy.visualization import simple_norm
from astropy.nddata import Cutout2D, NDData
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.table import Table, QTable
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.wcs.utils import pixel_to_skycoord
from astropy.coordinates import SkyCoord, match_coordinates_sky
from astropy.stats import sigma_clipped_stats
from photutils import CircularAperture, EPSFBuilder, find_peaks, CircularAnnulus
from photutils.detection import DAOStarFinder, IRAFStarFinder
from photutils.psf import DAOGroup, IntegratedGaussianPRF, extract_stars, IterativelySubtractedPSFPhotometry
from photutils.background import MMMBackground, MADStdBackgroundRMS
from photutils.centroids import centroid_2dg
from photutils import aperture_photometry
from ipywidgets import interact
import webbpsf
from webbpsf.utils import to_griddedpsfmodel
import pysynphot # PYSIN_CDBS must be defined in the user's environment (see note above)
/tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/lib/python3.7/site-packages/pysynphot/locations.py:345: UserWarning: Extinction files not found in /tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/grp/redcat/trds/extinction
warnings.warn('Extinction files not found in %s' % (extdir, ))
/tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/lib/python3.7/site-packages/pysynphot/refs.py:118: UserWarning: No graph or component tables found; functionality will be SEVERELY crippled. No files found for /tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/grp/redcat/trds/mtab/*_tmg.fits
'functionality will be SEVERELY crippled. ' + str(e))
/tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/lib/python3.7/site-packages/pysynphot/refs.py:125: UserWarning: No thermal tables found, no thermal calculations can be performed. No files found for /tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/grp/redcat/trds/mtab/*_tmt.fits
'no thermal calculations can be performed. ' + str(e))
**WARNING**: LOCAL JWST PRD VERSION PRDOPSSOC-034 DOESN'T MATCH THE CURRENT ONLINE VERSION PRDOPSSOC-039 Please consider updating pysiaf, e.g. pip install --upgrade pysiaf or conda update pysiaf
%matplotlib inline
from matplotlib import style, pyplot as plt
import matplotlib.patches as patches
import matplotlib.ticker as ticker
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 30
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 30
font1 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '12'}
font2 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '20'}
dict_images = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}
dict_filter_short = {}
dict_filter_long = {}
ff_short = []
det_short = []
det_long = []
ff_long = []
detlist_short = []
detlist_long = []
filtlist_short = []
filtlist_long = []
if not glob.glob('./*cal*fits'):
print("Downloading images")
boxlink_images_lev2 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level2.tar.gz'
boxfile_images_lev2 = './images_level2.tar.gz'
urllib.request.urlretrieve(boxlink_images_lev2, boxfile_images_lev2)
tar = tarfile.open(boxfile_images_lev2, 'r')
tar.extractall()
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))
else:
images_dir = './'
images = sorted(glob.glob(os.path.join(images_dir, "*cal.fits")))
for image in images:
im = fits.open(image)
f = im[0].header['FILTER']
d = im[0].header['DETECTOR']
if d == 'NRCBLONG':
d = 'NRCB5'
elif d == 'NRCALONG':
d = 'NRCA5'
else:
d = d
wv = np.float(f[1:3])
if wv > 24:
ff_long.append(f)
det_long.append(d)
else:
ff_short.append(f)
det_short.append(d)
detlist_short = sorted(list(dict.fromkeys(det_short)))
detlist_long = sorted(list(dict.fromkeys(det_long)))
unique_list_filters_short = []
unique_list_filters_long = []
for x in ff_short:
if x not in unique_list_filters_short:
dict_filter_short.setdefault(x, {})
for x in ff_long:
if x not in unique_list_filters_long:
dict_filter_long.setdefault(x, {})
for d_s in detlist_short:
dict_images[d_s] = dict_filter_short
for d_l in detlist_long:
dict_images[d_l] = dict_filter_long
filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))
if len(dict_images[d][f]) == 0:
dict_images[d][f] = {'images': [image]}
else:
dict_images[d][f]['images'].append(image)
print("Available Detectors for SW channel:", detlist_short)
print("Available Detectors for LW channel:", detlist_long)
print("Available SW Filters:", filtlist_short)
print("Available LW Filters:", filtlist_long)
Downloading images Available Detectors for SW channel: ['NRCB1'] Available Detectors for LW channel: ['NRCB5'] Available SW Filters: ['F070W', 'F115W', 'F200W'] Available LW Filters: ['F277W', 'F356W', 'F444W']
filters = ['F070W', 'F090W', 'F115W', 'F140M', 'F150W2', 'F150W', 'F162M', 'F164N', 'F182M',
'F187N', 'F200W', 'F210M', 'F212N', 'F250M', 'F277W', 'F300M', 'F322W2', 'F323N',
'F335M', 'F356W', 'F360M', 'F405N', 'F410M', 'F430M', 'F444W', 'F460M', 'F466N', 'F470N', 'F480M']
psf_fwhm = [0.987, 1.103, 1.298, 1.553, 1.628, 1.770, 1.801, 1.494, 1.990, 2.060, 2.141, 2.304, 2.341, 1.340,
1.444, 1.585, 1.547, 1.711, 1.760, 1.830, 1.901, 2.165, 2.179, 2.300, 2.302, 2.459, 2.507, 2.535, 2.574]
zp_modA = [25.7977, 25.9686, 25.8419, 24.8878, 27.0048, 25.6536, 24.6957, 22.3073, 24.8258, 22.1775, 25.3677, 24.3296,
22.1036, 22.7850, 23.5964, 24.8239, 23.6452, 25.3648, 20.8604, 23.5873, 24.3778, 23.4778, 20.5588,
23.2749, 22.3584, 23.9731, 21.9502, 20.0428, 19.8869, 21.9002]
zp_modB = [25.7568, 25.9771, 25.8041, 24.8738, 26.9821, 25.6279, 24.6767, 22.2903, 24.8042, 22.1499, 25.3391, 24.2909,
22.0574, 22.7596, 23.5011, 24.6792, 23.5769, 25.3455, 20.8631, 23.4885, 24.3883, 23.4555, 20.7007,
23.2763, 22.4677, 24.1562, 22.0422, 20.1430, 20.0173, 22.4086]
dict_utils = {filters[i]: {'psf fwhm': psf_fwhm[i], 'VegaMAG zp modA': zp_modA[i],
'VegaMAG zp modB': zp_modB[i]} for i in range(len(filters))}
dets_short = ['NRCB1'] # detector of interest in this example
filts_short = ['F115W', 'F200W'] # filters of interest in this example
# cell for display images using ipywidgets
def browse_images(images):
n = len(images)
def view_image(image):
det = 'NRCB1'
filt = 'F115W'
im = fits.open(dict_images[det][filt]['images'][image])
data_sb = im[1].data
norm = simple_norm(data_sb, 'sqrt', percent=99.)
plt.figure(figsize=(10, 10))
plt.title(filt)
plt.imshow(data_sb, norm=norm, cmap='Greys')
plt.show()
interact(view_image, image=(0, n - 1))
browse_images(dict_images['NRCB1']['F115W']['images'])
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(1, len(filts_short), i + 1)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(data_sb, norm=norm, cmap='Greys')
plt.tight_layout()
findfont: Font family ['helvetica'] not found. Falling back to DejaVu Sans.
dict_psfs_webbpsf = {}
for det in dets_short:
dict_psfs_webbpsf.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_psfs_webbpsf[det].setdefault(filt, {})
dict_psfs_webbpsf[det][filt]['psf model grid'] = None
dict_psfs_webbpsf[det][filt]['psf model single'] = None
def create_psf_model(fov=11, create_grid=False, num=9, save_psf=False, detsampled=False):
nrc = webbpsf.NIRCam()
nrc.detector = det
nrc.filter = filt
src = webbpsf.specFromSpectralType('G5V', catalog='phoenix')
if detsampled:
print("Creating a detector sampled PSF")
aa = 'detector sampled'
fov = 21
else:
print("Creating a oversampled PSF")
aa = 'oversampled'
fov = fov
print("Using a {field}".format(field=fov), "px fov")
if create_grid:
print("")
print("Creating a grid of PSF for filter {filt} and detector {det}".format(filt=filt, det=det))
print("")
num = num
if save_psf:
outname = "./PSF_%s_samp4_G5V_fov%d_npsfs%d.fits" % (filt, fov, num)
nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,
save=True, outfile=outname, use_detsampled_psf=detsampled)
else:
grid_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,
fov_pixels=fov, use_detsampled_psf=detsampled)
dict_psfs_webbpsf[det][filt]['psf model grid'] = grid_psf
else:
print("")
print("Creating a single PSF for filter {filt} and detector {det}".format(filt=filt, det=det))
print("")
num = 1
if save_psf:
outname = "./PSF_%s_samp4_G5V_fov%d_npsfs%d.fits" % (filt, fov, num)
nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,
save=True, outfile=outname, use_detsampled_psf=detsampled)
else:
single_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,
fov_pixels=fov, use_detsampled_psf=detsampled)
dict_psfs_webbpsf[det][filt]['psf model single'] = single_psf
return
for det in dets_short:
for filt in filts_short:
create_psf_model(fov=11, num=25, create_grid=False, save_psf=False, detsampled=False)
Creating a oversampled PSF
Using a 11 px fov
Creating a single PSF for filter F115W and detector NRCB1
Running instrument: NIRCam, filter: F115W
Running detector: NRCB1
Position 1/1: (1023, 1023) pixels
Creating a oversampled PSF
Using a 11 px fov
Creating a single PSF for filter F200W and detector NRCB1
Running instrument: NIRCam, filter: F200W
Running detector: NRCB1
Position 1/1: (1023, 1023) pixels
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
ax = plt.subplot(1, 2, i + 1)
norm_epsf = simple_norm(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], 'log', percent=99.)
ax.set_title(filt, fontsize=40)
ax.imshow(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], norm=norm_epsf)
ax.set_xlabel('X [px]', fontsize=30)
ax.set_ylabel('Y [px]', fontsize=30)
plt.tight_layout()
for det in dets_short:
for filt in filts_short:
create_psf_model(fov=11, num=25, create_grid=True, save_psf=False, detsampled=False)
Creating a oversampled PSF
Using a 11 px fov
Creating a grid of PSF for filter F115W and detector NRCB1
Running instrument: NIRCam, filter: F115W
Running detector: NRCB1
Position 1/25: (0, 0) pixels
Position 2/25: (0, 512) pixels
Position 3/25: (0, 1024) pixels
Position 4/25: (0, 1535) pixels
Position 5/25: (0, 2047) pixels
Position 6/25: (512, 0) pixels
Position 7/25: (512, 512) pixels
Position 8/25: (512, 1024) pixels
Position 9/25: (512, 1535) pixels
Position 10/25: (512, 2047) pixels
Position 11/25: (1024, 0) pixels
Position 12/25: (1024, 512) pixels
Position 13/25: (1024, 1024) pixels
Position 14/25: (1024, 1535) pixels
Position 15/25: (1024, 2047) pixels
Position 16/25: (1535, 0) pixels
Position 17/25: (1535, 512) pixels
Position 18/25: (1535, 1024) pixels
Position 19/25: (1535, 1535) pixels
Position 20/25: (1535, 2047) pixels
Position 21/25: (2047, 0) pixels
Position 22/25: (2047, 512) pixels
Position 23/25: (2047, 1024) pixels
Position 24/25: (2047, 1535) pixels
Position 25/25: (2047, 2047) pixels
Creating a oversampled PSF
Using a 11 px fov
Creating a grid of PSF for filter F200W and detector NRCB1
Running instrument: NIRCam, filter: F200W
Running detector: NRCB1
Position 1/25: (0, 0) pixels
Position 2/25: (0, 512) pixels
Position 3/25: (0, 1024) pixels
Position 4/25: (0, 1535) pixels
Position 5/25: (0, 2047) pixels
Position 6/25: (512, 0) pixels
Position 7/25: (512, 512) pixels
Position 8/25: (512, 1024) pixels
Position 9/25: (512, 1535) pixels
Position 10/25: (512, 2047) pixels
Position 11/25: (1024, 0) pixels
Position 12/25: (1024, 512) pixels
Position 13/25: (1024, 1024) pixels
Position 14/25: (1024, 1535) pixels
Position 15/25: (1024, 2047) pixels
Position 16/25: (1535, 0) pixels
Position 17/25: (1535, 512) pixels
Position 18/25: (1535, 1024) pixels
Position 19/25: (1535, 1535) pixels
Position 20/25: (1535, 2047) pixels
Position 21/25: (2047, 0) pixels
Position 22/25: (2047, 512) pixels
Position 23/25: (2047, 1024) pixels
Position 24/25: (2047, 1535) pixels
Position 25/25: (2047, 2047) pixels
webbpsf.gridded_library.display_psf_grid(dict_psfs_webbpsf[dets_short[0]][filts_short[0]]['psf model grid'],
zoom_in=False, figsize=(14, 14))
/tmp/nbcollection-ci/scanner-build-dir/jdat_notebooks/psf_photometry/lib/python3.7/site-packages/webbpsf/gridded_library.py:509: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it. axes[n-1-iy, ix].imshow(data[i], vmax=vmax, vmin=vmin, norm=norm)
dict_psfs_epsf = {}
for det in dets_short:
dict_psfs_epsf.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_psfs_epsf[det].setdefault(filt, {})
dict_psfs_epsf[det][filt]['table psf stars'] = {}
dict_psfs_epsf[det][filt]['epsf single'] = {}
dict_psfs_epsf[det][filt]['epsf grid'] = {}
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = None
dict_psfs_epsf[det][filt]['epsf single'][i + 1] = None
dict_psfs_epsf[det][filt]['epsf grid'][i + 1] = None
def find_stars_epsf(det='NRCA1', filt='F070W', dist_sel=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
print("Finding PSF stars on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px")
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th[j] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
psf_stars = daofind(data)
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = psf_stars
if dist_sel:
print("")
print("Calculating closest neigbhour distance")
d = []
daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
stars_tot = daofind_tot(data)
x_tot = stars_tot['xcentroid']
y_tot = stars_tot['ycentroid']
for xx, yy in zip(psf_stars['xcentroid'], psf_stars['ycentroid']):
sep = []
dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)
sep = np.sort(dist)[1:2][0]
d.append(sep)
psf_stars['min distance'] = d
mask_dist = (psf_stars['min distance'] > min_sep[j])
psf_stars = psf_stars[mask_dist]
dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = psf_stars
print("Minimum distance required:", min_sep[j], "px")
print("")
print("Number of isolated sources found in the image used to build ePSF for {f}:".format(f=filt), len(psf_stars))
print("-----------------------------------------------------")
print("")
else:
print("")
print("Number of sources used to build ePSF for {f}:".format(f=filt), len(psf_stars))
print("--------------------------------------------")
print("")
tic = time.perf_counter()
th = [700, 500] # threshold level for the two filters (length must match number of filters analyzed)
min_sep = [10, 10] # minimum separation acceptable for ePSF stars from closest neighbour
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
find_stars_epsf(det=det, filt=filt, dist_sel=False)
toc = time.perf_counter()
print("Elapsed Time for finding stars:", toc - tic)
Finding PSF stars on image 1 of filter F115W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 FWHM for the filter F115W: 1.298 px Number of sources used to build ePSF for F115W: 1306 -------------------------------------------- Finding PSF stars on image 2 of filter F115W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 FWHM for the filter F115W: 1.298 px Number of sources used to build ePSF for F115W: 1324 -------------------------------------------- Finding PSF stars on image 3 of filter F115W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 FWHM for the filter F115W: 1.298 px Number of sources used to build ePSF for F115W: 1307 -------------------------------------------- Finding PSF stars on image 4 of filter F115W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 FWHM for the filter F115W: 1.298 px Number of sources used to build ePSF for F115W: 1316 -------------------------------------------- Finding PSF stars on image 1 of filter F200W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 FWHM for the filter F200W: 2.141 px Number of sources used to build ePSF for F200W: 1276 -------------------------------------------- Finding PSF stars on image 2 of filter F200W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 FWHM for the filter F200W: 2.141 px Number of sources used to build ePSF for F200W: 1287 -------------------------------------------- Finding PSF stars on image 3 of filter F200W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 FWHM for the filter F200W: 2.141 px Number of sources used to build ePSF for F200W: 1291 -------------------------------------------- Finding PSF stars on image 4 of filter F200W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 FWHM for the filter F200W: 2.141 px Number of sources used to build ePSF for F200W: 1278 -------------------------------------------- Elapsed Time for finding stars: 32.79131748299983
def build_epsf(det='NRCA1', filt='F070W'):
mmm_bkg = MMMBackground()
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
data = data_sb / imh['PHOTMJSR']
hsize = (sizes[j] - 1) / 2
x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']
y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']
mask = ((x > hsize) & (x < (data.shape[1] - 1 - hsize)) & (y > hsize) & (y < (data.shape[0] - 1 - hsize)))
stars_tbl = Table()
stars_tbl['x'] = x[mask]
stars_tbl['y'] = y[mask]
bkg = mmm_bkg(data)
data_bkgsub = data.copy()
data_bkgsub -= bkg
nddata = NDData(data=data_bkgsub)
stars = extract_stars(nddata, stars_tbl, size=sizes[j])
print("Creating ePSF for image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
epsf_builder = EPSFBuilder(oversampling=oversample, maxiters=3, progress_bar=False)
epsf, fitted_stars = epsf_builder(stars)
dict_psfs_epsf[det][filt]['epsf single'][i + 1] = epsf
tic = time.perf_counter()
sizes = [11, 11] # size of the cutout (extract region) for each PSF star - must match number of filters analyzed
oversample = 4
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
build_epsf(det=det, filt=filt)
toc = time.perf_counter()
print("Time to build the Effective PSF:", toc - tic)
Creating ePSF for image 1 of filter F115W, detector NRCB1 Creating ePSF for image 2 of filter F115W, detector NRCB1 Creating ePSF for image 3 of filter F115W, detector NRCB1 Creating ePSF for image 4 of filter F115W, detector NRCB1 Creating ePSF for image 1 of filter F200W, detector NRCB1
WARNING: The star at (2024.6583712618988, 303.05613687467087) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1287.6938989443017, 184.04681924060748) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (2024.6583712618988, 303.05613687467087) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1385.8600398323995, 733.8972484228361) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (26.284959747891065, 926.5001627023972) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1942.2238557129556, 1446.3283575899338) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1356.7816503557121, 1694.7562979221789) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1913.8036243641827, 2036.2020641439606) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf]
Creating ePSF for image 2 of filter F200W, detector NRCB1
WARNING: The star at (1254.1534926597653, 37.344944331531984) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (2013.6414870311062, 312.81837167426096) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (31.360541297377488, 933.8276852413142) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (271.2395713411656, 1692.6433967443513) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1918.5987179055912, 2040.976602498569) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf]
Creating ePSF for image 3 of filter F200W, detector NRCB1
WARNING: The star at (811.4333523945055, 595.3465086660839) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1366.8241155624034, 231.93661176140287) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1366.715469929931, 232.74382880584332) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (2010.5398162755582, 315.48809568080753) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (378.35911922042163, 577.7662913926318) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (811.4333523945055, 595.3465086660839) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (970.297816824831, 816.1998334840716) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1212.6776286738234, 1638.0921925173056) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf]
Creating ePSF for image 4 of filter F200W, detector NRCB1
WARNING: The star at (975.8053091205506, 810.4505240833056) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1371.179005239543, 227.01540958771997) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1371.7935197529089, 227.0491003308669) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (2031.6502442013732, 304.50771304247894) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (2015.8337608261713, 309.5029288861399) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (380.0252804162881, 573.3385577797417) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (975.8053091205506, 810.4505240833056) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1252.4648604612446, 883.7542036883244) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf] WARNING: The star at (1217.2316895541658, 1632.3846435471307) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf]
Time to build the Effective PSF: 272.09734832799995
WARNING: The star at (1921.778778548655, 2037.860473681066) cannot be fit because its fitting region extends beyond the star cutout image. [photutils.psf.epsf]
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
ax = plt.subplot(1, 2, i + 1)
norm_epsf = simple_norm(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, 'log', percent=99.)
plt.title(filt, fontsize=30)
ax.imshow(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, norm=norm_epsf)
def count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40):
num_grid_calc = np.arange(2, grid_points + 1, 1)
num_grid_calc = num_grid_calc[::-1]
for num in num_grid_calc:
print("Calculating the number of PSF stars in a %d x %d grid:" % (num, num))
print("")
image = fits.open(dict_images[det][filt]['images'][i])
data_sb = image[1].data
points = np.int16((data_sb.shape[0] / num) / 2)
x_center = np.arange(points, 2 * points * (num), 2 * points)
y_center = np.arange(points, 2 * points * (num), 2 * points)
centers = np.array(np.meshgrid(x_center, y_center)).T.reshape(-1, 2)
for n, val in enumerate(centers):
x = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['xcentroid']
y = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['ycentroid']
flux = dict_psfs_epsf[det][filt]['table psf stars'][i + 1]['flux']
half_size = (size - 1) / 2
lim1 = val[0] - points + half_size
lim2 = val[0] + points - half_size
lim3 = val[1] - points + half_size
lim4 = val[1] + points - half_size
test = (x > lim1) & (x < lim2) & (y > lim3) & (y < lim4)
# if np.count_nonzero(test) < min_numpsf:
# raise ValueError("Not enough PSF stars in all the cells (> %d): Decrease your grid size or the minimum number of PSF stars in each cell or change parameters in the finder" %(min_numpsf))
if np.count_nonzero(test) < min_numpsf:
print("Center Coordinates of grid cell %d are (%d, %d) --- Not enough PSF stars in the cell (number of PSF stars < %d)" % (i + 1, val[0], val[1], min_numpsf))
else:
print("Center Coordinate of grid cell %d are (%d, %d) --- Number of PSF stars:" % (n + 1, val[0], val[1]), np.count_nonzero(test))
print("")
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
print("Analyzing image {number} of filter {f}, detector {d} ".format(number=i + 1, f=filt, d=det))
print("")
count_PSFstars_grid(grid_points=5, size=15, min_numpsf=40)
Analyzing image 1 of filter F115W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 72 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 46 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 51 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 48 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 51 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 55 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 51 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 45 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 40 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 57 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 47 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 50 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 55 Center Coordinates of grid cell 1 are (1020, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 15 are (1020, 1836) --- Number of PSF stars: 43 Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 50 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 44 Center Coordinates of grid cell 1 are (1428, 1020) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 47 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 50 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 45 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 50 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 54 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 44 Center Coordinate of grid cell 25 are (1836, 1836) --- Number of PSF stars: 40 Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 122 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 73 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 70 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 90 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 84 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 75 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 84 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 81 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 77 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 87 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 59 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 75 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 75 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 64 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 76 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 67 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 186 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 129 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 149 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 132 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 151 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 132 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 132 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 118 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 126 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 359 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 331 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 307 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 286 Analyzing image 2 of filter F115W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 74 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 49 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 50 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 45 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 48 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 59 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 53 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 43 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 62 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 50 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 49 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 55 Center Coordinate of grid cell 14 are (1020, 1428) --- Number of PSF stars: 42 Center Coordinate of grid cell 15 are (1020, 1836) --- Number of PSF stars: 42 Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 55 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 43 Center Coordinates of grid cell 2 are (1428, 1020) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 47 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 50 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 41 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 47 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 59 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 44 Center Coordinates of grid cell 2 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 138 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 79 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 75 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 82 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 84 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 73 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 84 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 86 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 80 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 82 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 61 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 69 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 72 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 62 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 70 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 73 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 201 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 134 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 150 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 139 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 143 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 133 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 126 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 123 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 125 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 378 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 333 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 303 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 280 Analyzing image 3 of filter F115W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 73 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 50 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 47 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 45 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 46 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 57 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 48 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 47 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 58 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 46 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 52 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 56 Center Coordinate of grid cell 14 are (1020, 1428) --- Number of PSF stars: 44 Center Coordinates of grid cell 3 are (1020, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 48 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 41 Center Coordinates of grid cell 3 are (1428, 1020) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 51 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 53 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 43 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 49 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 52 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 46 Center Coordinates of grid cell 3 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 128 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 77 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 69 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 81 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 81 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 73 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 87 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 79 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 75 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 82 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 65 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 74 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 71 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 65 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 75 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 69 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 191 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 127 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 147 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 135 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 140 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 133 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 131 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 120 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 131 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 366 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 320 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 297 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 292 Analyzing image 4 of filter F115W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 66 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 46 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 47 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 48 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 51 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 59 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 52 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 47 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 44 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 63 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 47 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 50 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 56 Center Coordinate of grid cell 14 are (1020, 1428) --- Number of PSF stars: 41 Center Coordinate of grid cell 15 are (1020, 1836) --- Number of PSF stars: 41 Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 54 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 45 Center Coordinate of grid cell 18 are (1428, 1020) --- Number of PSF stars: 40 Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 50 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 50 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 42 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 51 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 55 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 42 Center Coordinates of grid cell 4 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 121 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 76 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 69 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 90 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 86 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 76 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 84 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 83 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 79 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 89 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 63 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 71 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 72 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 67 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 71 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 70 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 188 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 131 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 155 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 138 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 145 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 132 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 134 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 121 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 124 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 366 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 333 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 312 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 281 Analyzing image 1 of filter F200W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 81 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 49 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 45 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 46 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 65 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 46 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 43 Center Coordinates of grid cell 1 are (612, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 57 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 43 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 47 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 55 Center Coordinates of grid cell 1 are (1020, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinates of grid cell 1 are (1020, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 46 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 41 Center Coordinate of grid cell 18 are (1428, 1020) --- Number of PSF stars: 41 Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 57 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 43 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 44 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 49 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 42 Center Coordinates of grid cell 1 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 149 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 79 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 65 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 81 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 75 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 68 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 80 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 77 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 71 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 89 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 57 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 71 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 72 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 56 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 70 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 69 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 211 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 121 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 142 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 127 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 145 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 123 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 123 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 107 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 128 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 376 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 308 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 291 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 277 Analyzing image 2 of filter F200W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 78 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 48 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 45 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 40 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 47 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 77 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 45 Center Coordinates of grid cell 2 are (612, 1020) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 43 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 60 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 47 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 48 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 56 Center Coordinates of grid cell 2 are (1020, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinates of grid cell 2 are (1020, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 49 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 42 Center Coordinate of grid cell 18 are (1428, 1020) --- Number of PSF stars: 42 Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 47 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 51 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 40 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 47 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 53 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 41 Center Coordinates of grid cell 2 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 155 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 77 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 63 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 80 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 80 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 69 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 81 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 80 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 77 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 91 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 57 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 67 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 68 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 58 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 66 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 67 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 216 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 119 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 142 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 132 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 147 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 124 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 124 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 110 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 121 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 385 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 309 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 301 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 263 Analyzing image 3 of filter F200W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 76 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 52 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 41 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 45 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 46 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 69 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 43 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 44 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 45 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 55 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 46 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 51 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 54 Center Coordinates of grid cell 3 are (1020, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinates of grid cell 3 are (1020, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 46 Center Coordinates of grid cell 3 are (1428, 612) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 18 are (1428, 1020) --- Number of PSF stars: 42 Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 55 Center Coordinate of grid cell 21 are (1836, 204) --- Number of PSF stars: 43 Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 45 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 51 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 43 Center Coordinates of grid cell 3 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 156 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 76 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 66 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 79 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 76 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 69 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 82 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 74 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 71 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 88 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 57 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 70 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 70 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 61 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 69 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 71 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 218 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 115 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 145 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 131 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 146 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 123 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 124 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 113 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 127 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 383 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 305 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 294 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 275 Analyzing image 4 of filter F200W, detector NRCB1 Calculating the number of PSF stars in a 5 x 5 grid: Center Coordinate of grid cell 1 are (204, 204) --- Number of PSF stars: 73 Center Coordinate of grid cell 2 are (204, 612) --- Number of PSF stars: 47 Center Coordinate of grid cell 3 are (204, 1020) --- Number of PSF stars: 42 Center Coordinate of grid cell 4 are (204, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 5 are (204, 1836) --- Number of PSF stars: 49 Center Coordinate of grid cell 6 are (612, 204) --- Number of PSF stars: 73 Center Coordinate of grid cell 7 are (612, 612) --- Number of PSF stars: 45 Center Coordinate of grid cell 8 are (612, 1020) --- Number of PSF stars: 42 Center Coordinate of grid cell 9 are (612, 1428) --- Number of PSF stars: 41 Center Coordinate of grid cell 10 are (612, 1836) --- Number of PSF stars: 57 Center Coordinate of grid cell 11 are (1020, 204) --- Number of PSF stars: 45 Center Coordinate of grid cell 12 are (1020, 612) --- Number of PSF stars: 45 Center Coordinate of grid cell 13 are (1020, 1020) --- Number of PSF stars: 54 Center Coordinates of grid cell 4 are (1020, 1428) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinates of grid cell 4 are (1020, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 16 are (1428, 204) --- Number of PSF stars: 52 Center Coordinate of grid cell 17 are (1428, 612) --- Number of PSF stars: 41 Center Coordinate of grid cell 18 are (1428, 1020) --- Number of PSF stars: 45 Center Coordinate of grid cell 19 are (1428, 1428) --- Number of PSF stars: 46 Center Coordinate of grid cell 20 are (1428, 1836) --- Number of PSF stars: 51 Center Coordinates of grid cell 4 are (1836, 204) --- Not enough PSF stars in the cell (number of PSF stars < 40) Center Coordinate of grid cell 22 are (1836, 612) --- Number of PSF stars: 49 Center Coordinate of grid cell 23 are (1836, 1020) --- Number of PSF stars: 55 Center Coordinate of grid cell 24 are (1836, 1428) --- Number of PSF stars: 42 Center Coordinates of grid cell 4 are (1836, 1836) --- Not enough PSF stars in the cell (number of PSF stars < 40) Calculating the number of PSF stars in a 4 x 4 grid: Center Coordinate of grid cell 1 are (256, 256) --- Number of PSF stars: 143 Center Coordinate of grid cell 2 are (256, 768) --- Number of PSF stars: 75 Center Coordinate of grid cell 3 are (256, 1280) --- Number of PSF stars: 64 Center Coordinate of grid cell 4 are (256, 1792) --- Number of PSF stars: 85 Center Coordinate of grid cell 5 are (768, 256) --- Number of PSF stars: 78 Center Coordinate of grid cell 6 are (768, 768) --- Number of PSF stars: 69 Center Coordinate of grid cell 7 are (768, 1280) --- Number of PSF stars: 80 Center Coordinate of grid cell 8 are (768, 1792) --- Number of PSF stars: 75 Center Coordinate of grid cell 9 are (1280, 256) --- Number of PSF stars: 76 Center Coordinate of grid cell 10 are (1280, 768) --- Number of PSF stars: 88 Center Coordinate of grid cell 11 are (1280, 1280) --- Number of PSF stars: 58 Center Coordinate of grid cell 12 are (1280, 1792) --- Number of PSF stars: 71 Center Coordinate of grid cell 13 are (1792, 256) --- Number of PSF stars: 68 Center Coordinate of grid cell 14 are (1792, 768) --- Number of PSF stars: 65 Center Coordinate of grid cell 15 are (1792, 1280) --- Number of PSF stars: 68 Center Coordinate of grid cell 16 are (1792, 1792) --- Number of PSF stars: 68 Calculating the number of PSF stars in a 3 x 3 grid: Center Coordinate of grid cell 1 are (341, 341) --- Number of PSF stars: 206 Center Coordinate of grid cell 2 are (341, 1023) --- Number of PSF stars: 117 Center Coordinate of grid cell 3 are (341, 1705) --- Number of PSF stars: 147 Center Coordinate of grid cell 4 are (1023, 341) --- Number of PSF stars: 129 Center Coordinate of grid cell 5 are (1023, 1023) --- Number of PSF stars: 145 Center Coordinate of grid cell 6 are (1023, 1705) --- Number of PSF stars: 122 Center Coordinate of grid cell 7 are (1705, 341) --- Number of PSF stars: 125 Center Coordinate of grid cell 8 are (1705, 1023) --- Number of PSF stars: 116 Center Coordinate of grid cell 9 are (1705, 1705) --- Number of PSF stars: 121 Calculating the number of PSF stars in a 2 x 2 grid: Center Coordinate of grid cell 1 are (512, 512) --- Number of PSF stars: 371 Center Coordinate of grid cell 2 are (512, 1536) --- Number of PSF stars: 310 Center Coordinate of grid cell 3 are (1536, 512) --- Number of PSF stars: 302 Center Coordinate of grid cell 4 are (1536, 1536) --- Number of PSF stars: 272
dict_phot = {}
for det in dets_short:
dict_phot.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_phot[det].setdefault(filt, {})
dict_phot[det][filt]['residual images'] = {}
dict_phot[det][filt]['output photometry tables'] = {}
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
dict_phot[det][filt]['residual images'][i + 1] = None
dict_phot[det][filt]['output photometry tables'][i + 1] = None
def psf_phot(det='NRCA1', filt='F070W', th=2000, psf='grid_webbpsf', save_residuals=False, save_output=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
fitter = LevMarLSQFitter()
im = fits.open(dict_images[det][filt]['images'][i])
imh = im[1].header
data_sb = im[1].data
d = im[0].header['DETECTOR']
prim_dith_pos = im[0].header['PATT_NUM']
prim_dith_num = im[0].header['NUMDTHPT']
subpx_dith_pos = im[0].header['SUBPXNUM']
subpx_dith_num = im[0].header['SUBPXPNS']
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
print("Applying conversion to the data")
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf)
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
daogroup = DAOGroup(5.0 * sigma_psf)
# grid PSF
if psf == 'grid_webbpsf':
print("Using as PSF model WebbPSF PSFs grid")
psf_model = dict_psfs_webbpsf[det][filt]['psf model grid'].copy()
# single psf:
if psf == 'single_webbpsf':
print("Using as PSF model WebbPSF single PSF")
psf_model = dict_psfs_webbpsf[det][filt]['psf model single'].copy()
# epsf:
if psf == 'single_epsf':
print("Using as PSF model single ePSF")
psf_model = dict_psfs_epsf[det][filt]['epsf single'][i + 1].copy()
print("Performing the photometry on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
tic = time.perf_counter()
phot = IterativelySubtractedPSFPhotometry(finder=daofind, group_maker=daogroup,
bkg_estimator=mmm_bkg, psf_model=psf_model,
fitter=LevMarLSQFitter(),
niters=2, fitshape=(11, 11), aperture_radius=ap_radius[j])
result = phot(data)
toc = time.perf_counter()
print("Time needed to perform photometry on image {number}:".format(number=i + 1), "%.2f" % ((toc - tic) / 3600), "hours")
print("Number of sources detected in image {number} for filter {f}:".format(number=i + 1, f=filt), len(result))
residual_image = phot.get_residual_image()
dict_phot[det][filt]['residual images'][i + 1] = residual_image
dict_phot[det][filt]['output photometry tables'][i + 1] = result
# save the residual images as fits file:
if save_residuals:
hdu = fits.PrimaryHDU(residual_image)
hdul = fits.HDUList([hdu])
residual_outname = 'residual_%s_%s_webbPSF_gridPSF_%dof%d_%dof%d.fits' % (d, filt, prim_dith_pos, prim_dith_num, subpx_dith_pos, subpx_dith_num)
dir_output_phot = './'
hdul.writeto(os.path.join(dir_output_phot, residual_outname))
outname = 'phot_%s_%s_webbPSF_gridPSF_level2_%dof%d_%dof%d.pkl' % (d, filt, prim_dith_pos, prim_dith_num, subpx_dith_pos, subpx_dith_num)
# save the output photometry Tables
if save_output:
tab = result.to_pandas()
tab.to_pickle(os.path.join(dir_output_phot, outname))
tic_tot = time.perf_counter()
ap_radius = [3.0, 3.5] # must match the number of filters analyzed
if glob.glob('./*residual*.fits'):
print("Deleting Residual images from directory")
files = glob.glob('./residual*.fits')
for file in files:
os.remove(file)
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images[det][filt]['images']), 1):
psf_phot(det=det, filt=filt, th=2000, psf='grid_webbpsf', save_residuals=True, save_output=False)
toc_tot = time.perf_counter()
print("Time elapsed to perform the photometry of the {number} images:".format(number=(len(filts_short) * len(dict_images[det][filt]['images']))), "%.2f" % ((toc_tot - tic_tot) / 3600), "hours")
Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 Applying conversion to the data FWHM for the filter F115W: 1.298 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 1 of filter F115W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 1: 0.01 hours Number of sources detected in image 1 for filter F115W: 557 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 Applying conversion to the data FWHM for the filter F115W: 1.298 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 2 of filter F115W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 2: 0.01 hours Number of sources detected in image 2 for filter F115W: 567 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 Applying conversion to the data FWHM for the filter F115W: 1.298 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 3 of filter F115W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 3: 0.02 hours Number of sources detected in image 3 for filter F115W: 587 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 Applying conversion to the data FWHM for the filter F115W: 1.298 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 4 of filter F115W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 4: 0.01 hours Number of sources detected in image 4 for filter F115W: 566 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 Applying conversion to the data FWHM for the filter F200W: 2.141 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 1 of filter F200W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 1: 0.01 hours Number of sources detected in image 1 for filter F200W: 402 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 Applying conversion to the data FWHM for the filter F200W: 2.141 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 2 of filter F200W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 2: 0.01 hours Number of sources detected in image 2 for filter F200W: 428 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 Applying conversion to the data FWHM for the filter F200W: 2.141 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 3 of filter F200W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 3: 0.01 hours Number of sources detected in image 3 for filter F200W: 422 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 Applying conversion to the data FWHM for the filter F200W: 2.141 Using as PSF model WebbPSF PSFs grid Performing the photometry on image 4 of filter F200W, detector NRCB1
WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting] WARNING: The fit may be unsuccessful; check fit_info['message'] for more information. [astropy.modeling.fitting]
Time needed to perform photometry on image 4: 0.01 hours Number of sources detected in image 4 for filter F200W: 428 Time elapsed to perform the photometry of the 8 images: 0.11 hours
dict_phot['NRCB1']['F115W']['output photometry tables'][1]
| x_0 | x_fit | y_0 | y_fit | flux_0 | flux_fit | id | group_id | flux_unc | x_0_unc | y_0_unc | iter_detected |
|---|---|---|---|---|---|---|---|---|---|---|---|
| float64 | float64 | float64 | float64 | float64 | float64 | int64 | int64 | float64 | float64 | float64 | int64 |
| 319.4246272573411 | 319.6243829607715 | 16.698703468391624 | 16.779812043448285 | 501.7585915115253 | 817.362838092945 | 1 | 1 | 8.990221716253659 | 0.009307183796375833 | 0.010931316746915665 | 1 |
| 239.07687209767755 | 239.0617950028212 | 21.59248084751962 | 21.389574086731667 | 343.62489297091554 | 557.6927464772607 | 2 | 2 | 6.126577686641626 | 0.011929795802084938 | 0.009337876774809465 | 1 |
| 486.6367667247977 | 486.39907839623015 | 25.298505419946128 | 25.581923420634563 | 604.8895459014808 | 986.9494366579039 | 3 | 3 | 10.83111053252194 | 0.00899513797852723 | 0.009309945475638176 | 1 |
| 1498.583912593049 | 1498.3732082503427 | 25.6631229860481 | 25.767128766572036 | 378.6391370981348 | 615.3356760509652 | 4 | 4 | 6.742970761608649 | 0.00907145032403124 | 0.010986404963736896 | 1 |
| 1772.225384572788 | 1772.1571288023129 | 25.870632629745227 | 25.91433319889012 | 309.433462977977 | 502.15276048146166 | 5 | 5 | 5.417956435341972 | 0.010304378070937535 | 0.012247849065957368 | 1 |
| 1249.3889857374545 | 1249.7823237777382 | 33.16566796574103 | 33.33673282204046 | 5405.539907386959 | 4381.280446601617 | 6 | 6 | 391.0343000665844 | 0.09184877509422591 | 0.07949213846486462 | 1 |
| 1969.2495645233275 | 1969.162796362333 | 33.29229202776429 | 33.58300446690942 | 867.9795565502349 | 1414.6972378229914 | 7 | 7 | 15.144823308239118 | 0.010144646990381765 | 0.00912842297131713 | 1 |
| 1815.8771749967136 | 1815.907337301414 | 39.23846441749301 | 39.28437544313809 | 2717.042432560357 | 3547.9737199442006 | 8 | 8 | 100.09189430852692 | 0.03160034295248331 | 0.025917024937055357 | 1 |
| 524.6487415112808 | 524.406998298359 | 47.02432602926747 | 47.030845271298716 | 472.5076338924839 | 766.2118980024395 | 9 | 9 | 8.203052068604297 | 0.008814536490983032 | 0.012107564445447182 | 1 |
| 283.71266677261855 | 283.7222137848495 | 51.77102831013282 | 51.74036806782724 | 1290.3598304109435 | 1732.2634939523075 | 10 | 10 | 42.74888967212685 | 0.02252993269319067 | 0.023640700373578547 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 313.86603305995726 | 313.19682506135877 | 1514.5427673094175 | 1514.8563278296463 | 342.1082956933899 | 603.8680259453085 | 33 | 25 | 75.23217712488908 | 0.12239568801783855 | 0.14170849173789118 | 2 |
| 1594.8956462726592 | 1595.1250833802517 | 1548.8100481748231 | 1548.9063013229852 | 472.58088565049104 | 629.8278438885151 | 34 | 26 | 62.19146680758922 | 0.09713292118206537 | 0.11201900873414859 | 2 |
| 1662.1578859919007 | 1663.2265437718502 | 1707.1181529961095 | 1707.2296299090353 | -1152.4314388705366 | -801.7816349532864 | 35 | 27 | 220.42205001600678 | 0.24848655824329313 | 0.246330399739303 | 2 |
| 887.2699332861704 | 887.3010483357664 | 1771.8272770710769 | 1771.814971784264 | 414.8638966898997 | 689.2459645066103 | 36 | 28 | 50.941044736391476 | 0.06538070260502926 | 0.08249922169550275 | 2 |
| 632.1139545087674 | 632.2413361589065 | 1832.131903961262 | 1833.319054508187 | -1209.7248563610865 | -869.4322728076429 | 37 | 29 | 192.42715994826693 | 0.20475413236066448 | 0.1883377953147412 | 2 |
| 1953.5607278538305 | 1953.7049301558998 | 1878.2066905345011 | 1878.613283434624 | 373.2331683449073 | 669.9881881867456 | 38 | 30 | 81.11772983168039 | 0.1251978399831555 | 0.11795975423024868 | 2 |
| 81.00891600674758 | 80.81134851132327 | 1946.952614338015 | 1948.278489889208 | -1219.1934508690745 | -919.9909476407366 | 39 | 31 | 184.34694394562416 | 0.19686041888868502 | 0.17335198358656354 | 2 |
| 220.02187625848927 | 220.2836683601479 | 1946.9740884823827 | 1948.2627148391628 | -904.5799652848325 | -679.9729307831898 | 40 | 32 | 137.33832585373585 | 0.1838440155627371 | 0.1755483860355323 | 2 |
| 904.6863753881605 | 904.1933825155567 | 1981.665271054273 | 1981.8892865120495 | 399.6678295172721 | 615.9851749936068 | 41 | 33 | 67.32396610431962 | 0.10427308333664896 | 0.12267888876663632 | 2 |
| 1914.8493639376886 | 1914.3952207936882 | 2037.9655624575653 | 2037.2523008520304 | 785.576598102754 | 1430.7031270837836 | 42 | 34 | 294.07734075602554 | 0.1787639967246127 | 0.18490724530635988 | 2 |
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(2, len(filts_short), i + 1)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(data_sb, norm=norm, cmap='Greys')
for det in dets_short:
for i, filt in enumerate(filts_short):
res = dict_phot[det][filt]['residual images'][1]
ax = plt.subplot(2, len(filts_short), i + 3)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
ax.imshow(res, norm=norm, cmap='Greys')
plt.tight_layout()
if not glob.glob('./*phot*gridPSF*.pkl'):
print("Downloading Photometry Output")
boxlink_cat_f115w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F115W.tar.gz'
boxfile_cat_f115w = './phot_cat_F115W.tar.gz'
urllib.request.urlretrieve(boxlink_cat_f115w, boxfile_cat_f115w)
tar = tarfile.open(boxfile_cat_f115w, 'r')
tar.extractall()
boxlink_cat_f200w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F200W.tar.gz'
boxfile_cat_f200w = './phot_cat_F200W.tar.gz'
urllib.request.urlretrieve(boxlink_cat_f200w, boxfile_cat_f200w)
tar = tarfile.open(boxfile_cat_f200w, 'r')
tar.extractall()
cat_dir = './'
phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))
else:
cat_dir = './'
phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))
phots_pkl_f200w = sorted(glob.glob(os.path.join(cat_dir, '*F200W*gridPSF*.pkl')))
results_f115w = []
results_f200w = []
for phot_pkl_f115w, phot_pkl_f200w in zip(phots_pkl_f115w, phots_pkl_f200w):
ph_f115w = pd.read_pickle(phot_pkl_f115w)
ph_f200w = pd.read_pickle(phot_pkl_f200w)
result_f115w = QTable.from_pandas(ph_f115w)
result_f200w = QTable.from_pandas(ph_f200w)
results_f115w.append(result_f115w)
results_f200w.append(result_f200w)
Downloading Photometry Output
images_f115w = []
images_f200w = []
for i in np.arange(0, len(dict_images['NRCB1']['F115W']['images']), 1):
image_f115w = ImageModel(dict_images['NRCB1']['F115W']['images'][i])
images_f115w.append(image_f115w)
for i in np.arange(0, len(dict_images['NRCB1']['F200W']['images']), 1):
image_f200w = ImageModel(dict_images['NRCB1']['F200W']['images'][i])
images_f200w.append(image_f200w)
results_clean_f115w = []
results_clean_f200w = []
for i in np.arange(0, len(images_f115w), 1):
mask_f115w = ((results_f115w[i]['x_fit'] > 0) & (results_f115w[i]['x_fit'] < 2048) &
(results_f115w[i]['y_fit'] > 0) & (results_f115w[i]['y_fit'] < 2048) &
(results_f115w[i]['flux_fit'] > 0))
result_clean_f115w = results_f115w[i][mask_f115w]
ra_f115w, dec_f115w = images_f115w[i].meta.wcs(result_clean_f115w['x_fit'], result_clean_f115w['y_fit'])
radec_f115w = SkyCoord(ra_f115w, dec_f115w, unit='deg')
result_clean_f115w['radec'] = radec_f115w
results_clean_f115w.append(result_clean_f115w)
mask_f200w = ((results_f200w[i]['x_fit'] > 0) & (results_f200w[i]['x_fit'] < 2048) &
(results_f200w[i]['y_fit'] > 0) & (results_f200w[i]['y_fit'] < 2048) &
(results_f200w[i]['flux_fit'] > 0))
result_clean_f200w = results_f200w[i][mask_f200w]
ra_f200w, dec_f200w = images_f200w[i].meta.wcs(result_clean_f200w['x_fit'], result_clean_f200w['y_fit'])
radec_f200w = SkyCoord(ra_f200w, dec_f200w, unit='deg')
result_clean_f200w['radec'] = radec_f200w
results_clean_f200w.append(result_clean_f200w)
max_sep = 0.015 * u.arcsec
matches_phot_single = []
filt1 = 'F115W'
filt2 = 'F200W'
for res1, res2 in zip(results_clean_f115w, results_clean_f200w):
idx, d2d, _ = match_coordinates_sky(res1['radec'], res2['radec'])
sep_constraint = d2d < max_sep
match_phot_single = Table()
x_0_f115w = res1['x_0'][sep_constraint]
y_0_f115w = res1['y_0'][sep_constraint]
x_fit_f115w = res1['x_fit'][sep_constraint]
y_fit_f115w = res1['y_fit'][sep_constraint]
radec_f115w = res1['radec'][sep_constraint]
mag_f115w = (-2.5 * np.log10(res1['flux_fit']))[sep_constraint]
emag_f115w = (1.086 * (res1['flux_unc'] / res1['flux_fit']))[sep_constraint]
x_0_f200w = res2['x_0'][idx[sep_constraint]]
y_0_f200w = res2['y_0'][idx[sep_constraint]]
x_fit_f200w = res2['x_fit'][idx[sep_constraint]]
y_fit_f200w = res2['y_fit'][idx[sep_constraint]]
radec_f200w = res2['radec'][idx][sep_constraint]
mag_f200w = (-2.5 * np.log10(res2['flux_fit']))[idx[sep_constraint]]
emag_f200w = (1.086 * (res2['flux_unc'] / res2['flux_fit']))[idx[sep_constraint]]
match_phot_single['x_0_' + filt1] = x_0_f115w
match_phot_single['y_0_' + filt1] = y_0_f115w
match_phot_single['x_fit_' + filt1] = x_fit_f115w
match_phot_single['y_fit_' + filt1] = y_fit_f115w
match_phot_single['radec_' + filt1] = radec_f115w
match_phot_single['mag_' + filt1] = mag_f115w
match_phot_single['emag_' + filt1] = emag_f115w
match_phot_single['x_0_' + filt2] = x_0_f200w
match_phot_single['y_0_' + filt2] = y_0_f200w
match_phot_single['x_fit_' + filt2] = x_fit_f200w
match_phot_single['y_fit_' + filt2] = y_fit_f200w
match_phot_single['radec_' + filt2] = radec_f200w
match_phot_single['mag_' + filt2] = mag_f200w
match_phot_single['emag_' + filt2] = emag_f200w
matches_phot_single.append(match_phot_single)
plt.figure(figsize=(12, 16))
plt.clf()
for i in np.arange(0, len(matches_phot_single), 1):
ax = plt.subplot(2, 2, i + 1)
j = str(i + 1)
xlim0 = -0.5
xlim1 = 0.8
ylim0 = -1
ylim1 = -9
ax.set_xlim(xlim0, xlim1)
ax.set_ylim(ylim0, ylim1)
ax.xaxis.set_major_locator(ticker.AutoLocator())
ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax.yaxis.set_major_locator(ticker.AutoLocator())
ax.yaxis.set_minor_locator(ticker.AutoMinorLocator())
f115w_single = matches_phot_single[i]['mag_' + filt1]
f200w_single = matches_phot_single[i]['mag_' + filt2]
ax.scatter(f115w_single - f200w_single, f115w_single, s=1, color='k')
ax.set_xlabel(filt1 + '-' + filt2, fontdict=font2)
ax.set_ylabel(filt1, fontdict=font2)
ax.text(xlim0 + 0.1, -8.65, "Image %s" % j, fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(12, 6))
ax1 = plt.subplot(1, 2, 1)
xlim0 = -1
xlim1 = 1
ylim0 = -1
ylim1 = 1
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
x_find_f115w = results_clean_f115w[0]['x_0']
y_find_f115w = results_clean_f115w[0]['y_0']
x_psf_f115w = results_clean_f115w[0]['x_fit']
y_psf_f115w = results_clean_f115w[0]['y_fit']
delta_x_f115w = x_find_f115w - x_psf_f115w
delta_y_f115w = y_find_f115w - y_psf_f115w
_, d_x_f115w, sigma_d_x_f115w = sigma_clipped_stats(delta_x_f115w)
_, d_y_f115w, sigma_d_y_f115w = sigma_clipped_stats(delta_y_f115w)
ax1.scatter(delta_x_f115w, delta_y_f115w, s=1, color='gray')
ax1.set_xlabel('$\Delta$ X (px)', fontdict=font2)
ax1.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax1.set_title(filt1, fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 0.15, ' $\Delta$ X = %5.3f $\pm$ %5.3f' % (d_x_f115w, sigma_d_x_f115w),
color='k', fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 0.30, ' $\Delta$ Y = %5.3f $\pm$ %5.3f' % (d_y_f115w, sigma_d_y_f115w),
color='k', fontdict=font2)
ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2 = plt.subplot(1, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
x_find_f200w = results_clean_f200w[0]['x_0']
y_find_f200w = results_clean_f200w[0]['y_0']
x_psf_f200w = results_clean_f200w[0]['x_fit']
y_psf_f200w = results_clean_f200w[0]['y_fit']
delta_x_f200w = x_find_f200w - x_psf_f200w
delta_y_f200w = y_find_f200w - y_psf_f200w
_, d_x_f200w, sigma_d_x_f200w = sigma_clipped_stats(delta_x_f200w)
_, d_y_f200w, sigma_d_y_f200w = sigma_clipped_stats(delta_y_f200w)
ax2.scatter(delta_x_f200w, delta_y_f200w, s=1, color='gray')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, ' $\Delta$ X = %5.3f $\pm$ %5.3f' % (d_x_f200w, sigma_d_x_f200w),
color='k', fontdict=font2)
ax2.text(xlim0 + 0.05, ylim1 - 0.30, ' $\Delta$ Y = %5.3f $\pm$ %5.3f' % (d_y_f200w, sigma_d_y_f200w),
color='k', fontdict=font2)
ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.set_xlabel('$\Delta$ X (px)', fontdict=font2)
ax2.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax2.set_title(filt2, fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(12, 8))
ax1 = plt.subplot(2, 2, 1)
xlim0 = -9
xlim1 = -1
ylim0 = -1
ylim1 = 1
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
mag_inst_f115w = -2.5 * np.log10(results_clean_f115w[0]['flux_fit'])
ax1.scatter(mag_inst_f115w, delta_x_f115w, s=1, color='gray')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax1.set_xlabel(filt1 + '_inst', fontdict=font2)
ax1.set_ylabel('$\Delta$ X (px)', fontdict=font2)
ax2 = plt.subplot(2, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(mag_inst_f115w, delta_y_f115w, s=1, color='gray')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.set_xlabel(filt1 + '_inst', fontdict=font2)
ax2.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
ax3 = plt.subplot(2, 2, 3)
ax3.set_xlim(xlim0, xlim1)
ax3.set_ylim(ylim0, ylim1)
ax3.xaxis.set_major_locator(ticker.AutoLocator())
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.yaxis.set_major_locator(ticker.AutoLocator())
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())
mag_inst_f200w = -2.5 * np.log10(results_clean_f200w[0]['flux_fit'])
ax3.scatter(mag_inst_f200w, delta_x_f200w, s=1, color='gray')
ax3.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax3.set_xlabel(filt2 + '_inst', fontdict=font2)
ax3.set_ylabel('$\Delta$ X (px)', fontdict=font2)
ax4 = plt.subplot(2, 2, 4)
ax4.set_xlim(xlim0, xlim1)
ax4.set_ylim(ylim0, ylim1)
ax4.xaxis.set_major_locator(ticker.AutoLocator())
ax4.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax4.yaxis.set_major_locator(ticker.AutoLocator())
ax4.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax4.scatter(mag_inst_f200w, delta_y_f200w, s=1, color='gray')
ax4.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax4.set_xlabel(filt2 + '_inst', fontdict=font2)
ax4.set_ylabel('$\Delta$ Y (px)', fontdict=font2)
plt.tight_layout()
def crossmatch_filter(table=None):
num = 0
num_cat = np.char.mod('%d', np.arange(1, len(table) + 1, 1))
idx_12, d2d_12, _ = match_coordinates_sky(table[num]['radec'], table[num + 1]['radec'])
sep_constraint_12 = d2d_12 < max_sep
matches_12 = Table()
matches_12['radec_' + num_cat[num]] = table[num]['radec'][sep_constraint_12]
matches_12['mag_' + num_cat[num]] = (-2.5 * np.log10(table[num]['flux_fit']))[sep_constraint_12]
matches_12['emag_' + num_cat[num]] = (1.086 * (table[num]['flux_unc'] /
table[num]['flux_fit']))[sep_constraint_12]
matches_12['radec_' + num_cat[num + 1]] = table[num + 1]['radec'][idx_12[sep_constraint_12]]
matches_12['mag_' + num_cat[num + 1]] = (-2.5 * np.log10(table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]
matches_12['emag_' + num_cat[num + 1]] = (1.086 * (table[num + 1]['flux_unc'] /
table[num + 1]['flux_fit']))[idx_12[sep_constraint_12]]
idx_123, d2d_123, _ = match_coordinates_sky(matches_12['radec_' + num_cat[num]], table[num + 2]['radec'])
sep_constraint_123 = d2d_123 < max_sep
matches_123 = Table()
matches_123['radec_' + num_cat[num]] = matches_12['radec_' + num_cat[num]][sep_constraint_123]
matches_123['mag_' + num_cat[num]] = matches_12['mag_' + num_cat[num]][sep_constraint_123]
matches_123['emag_' + num_cat[num]] = matches_12['emag_' + num_cat[num]][sep_constraint_123]
matches_123['radec_' + num_cat[num + 1]] = matches_12['radec_' + num_cat[num + 1]][sep_constraint_123]
matches_123['mag_' + num_cat[num + 1]] = matches_12['mag_' + num_cat[num + 1]][sep_constraint_123]
matches_123['emag_' + num_cat[num + 1]] = matches_12['emag_' + num_cat[num + 1]][sep_constraint_123]
matches_123['radec_' + num_cat[num + 2]] = table[num + 2]['radec'][idx_123[sep_constraint_123]]
matches_123['mag_' + num_cat[num + 2]] = (-2.5 * np.log10(table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]
matches_123['emag_' + num_cat[num + 2]] = (1.086 * (table[num + 2]['flux_unc'] /
table[num + 2]['flux_fit']))[idx_123[sep_constraint_123]]
idx_1234, d2d_1234, _ = match_coordinates_sky(matches_123['radec_' + num_cat[num]], table[num + 3]['radec'])
sep_constraint_1234 = d2d_1234 < max_sep
matches_1234 = Table()
matches_1234['radec_' + num_cat[num]] = matches_123['radec_' + num_cat[num]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num]] = matches_123['mag_' + num_cat[num]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num]] = matches_123['emag_' + num_cat[num]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 1]] = matches_123['radec_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num + 1]] = matches_123['mag_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num + 1]] = matches_123['emag_' + num_cat[num + 1]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 2]] = matches_123['radec_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['mag_' + num_cat[num + 2]] = matches_123['mag_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['emag_' + num_cat[num + 2]] = matches_123['emag_' + num_cat[num + 2]][sep_constraint_1234]
matches_1234['radec_' + num_cat[num + 3]] = table[num + 3]['radec'][idx_1234[sep_constraint_1234]]
matches_1234['mag_' + num_cat[num + 3]] = (-2.5 * np.log10(table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]
matches_1234['emag_' + num_cat[num + 3]] = (1.086 * (table[num + 3]['flux_unc'] /
table[num + 3]['flux_fit']))[idx_1234[sep_constraint_1234]]
matches_1234
return matches_1234
matches_f115w = crossmatch_filter(table=results_clean_f115w)
matches_f200w = crossmatch_filter(table=results_clean_f200w)
df_f115w = matches_f115w.to_pandas()
df_f200w = matches_f200w.to_pandas()
df_f115w['RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)
df_f115w['e_RA_' + filt1] = df_f115w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)
df_f115w['Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)
df_f115w['e_Dec_' + filt1] = df_f115w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)
df_f115w[filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)
df_f115w['e' + filt1 + '_inst'] = df_f115w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)
df_f200w['RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].mean(axis=1)
df_f200w['e_RA_' + filt2] = df_f200w[['radec_1.ra', 'radec_2.ra', 'radec_3.ra', 'radec_4.ra']].std(axis=1)
df_f200w['Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].mean(axis=1)
df_f200w['e_Dec_' + filt2] = df_f200w[['radec_1.dec', 'radec_2.dec', 'radec_3.dec', 'radec_4.dec']].std(axis=1)
df_f200w[filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].mean(axis=1)
df_f200w['e' + filt2 + '_inst'] = df_f200w[['mag_1', 'mag_2', 'mag_3', 'mag_4']].std(axis=1)
plt.figure(figsize=(12, 14))
plt.clf()
ax1 = plt.subplot(1, 2, 1)
ax1.set_xlabel(filt1 + '_inst -' + filt2 + '_inst', fontdict=font2)
ax1.set_ylabel(filt1 + '_inst', fontdict=font2)
xlim0 = -0.5
xlim1 = 0.8
ylim0 = -1.5
ylim1 = -9
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
radec_f115w_inst = SkyCoord(df_f115w['RA_' + filt1], df_f115w['Dec_' + filt1], unit='deg')
radec_f200w_inst = SkyCoord(df_f200w['RA_' + filt2], df_f200w['Dec_' + filt2], unit='deg')
idx_inst, d2d_inst, _ = match_coordinates_sky(radec_f115w_inst, radec_f200w_inst)
sep_constraint_inst = d2d_inst < max_sep
f115w_inst = np.array(df_f115w[filt1 + '_inst'][sep_constraint_inst])
ef115w_inst = np.array(df_f115w['e' + filt1 + '_inst'][sep_constraint_inst])
radec_f115w = radec_f115w_inst[sep_constraint_inst]
f200w_inst = np.array(df_f200w[filt2 + '_inst'][idx_inst[sep_constraint_inst]])
ef200w_inst = np.array(df_f200w['e' + filt2 + '_inst'][idx_inst[sep_constraint_inst]])
radec_f200w = radec_f200w_inst[idx_inst[sep_constraint_inst]]
ax1.scatter(f115w_inst - f200w_inst, f115w_inst, s=1, color='k')
ax2 = plt.subplot(2, 2, 2)
ax2.set_xlabel(filt1 + '_inst', fontdict=font2)
ax2.set_ylabel('$\sigma$' + filt1, fontdict=font2)
xlim0 = -9
xlim1 = -1.5
ylim0 = -0.01
ylim1 = 1
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(df_f115w[filt1 + '_inst'], df_f115w['e' + filt1 + '_inst'], s=1, color='k')
ax3 = plt.subplot(2, 2, 4)
ax3.set_xlabel(filt2 + '_inst', fontdict=font2)
ax3.set_ylabel('$\sigma$' + filt2, fontdict=font2)
ax3.set_xlim(xlim0, xlim1)
ax3.set_ylim(ylim0, ylim1)
ax3.xaxis.set_major_locator(ticker.AutoLocator())
ax3.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.yaxis.set_major_locator(ticker.AutoLocator())
ax3.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax3.scatter(df_f200w[filt2 + '_inst'], df_f200w['e' + filt2 + '_inst'], s=1, color='k')
plt.tight_layout()
dict_images_combined = {'NRCA1': {}, 'NRCA2': {}, 'NRCA3': {}, 'NRCA4': {}, 'NRCA5': {},
'NRCB1': {}, 'NRCB2': {}, 'NRCB3': {}, 'NRCB4': {}, 'NRCB5': {}}
dict_filter_short = {}
dict_filter_long = {}
ff_short = []
det_short = []
det_long = []
ff_long = []
detlist_short = []
detlist_long = []
filtlist_short = []
filtlist_long = []
if not glob.glob('./*combined*fits'):
print("Downloading images")
boxlink_images_lev3 = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/images_level3.tar.gz'
boxfile_images_lev3 = './images_level3.tar.gz'
urllib.request.urlretrieve(boxlink_images_lev3, boxfile_images_lev3)
tar = tarfile.open(boxfile_images_lev3, 'r')
tar.extractall()
images_dir = './'
files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
else:
images_dir = './'
files_singles = sorted(glob.glob(os.path.join(images_dir, "*combined*fits")))
for file in files_singles:
im = fits.open(file)
f = im[0].header['FILTER']
d = im[0].header['DETECTOR']
if d == 'NRCBLONG':
d = 'NRCB5'
elif d == 'NRCALONG':
d = 'NRCA5'
else:
d = d
wv = np.float(f[1:3])
if wv > 24:
ff_long.append(f)
det_long.append(d)
else:
ff_short.append(f)
det_short.append(d)
detlist_short = sorted(list(dict.fromkeys(det_short)))
detlist_long = sorted(list(dict.fromkeys(det_long)))
unique_list_filters_short = []
unique_list_filters_long = []
for x in ff_short:
if x not in unique_list_filters_short:
dict_filter_short.setdefault(x, {})
for x in ff_long:
if x not in unique_list_filters_long:
dict_filter_long.setdefault(x, {})
for d_s in detlist_short:
dict_images_combined[d_s] = dict_filter_short
for d_l in detlist_long:
dict_images_combined[d_l] = dict_filter_long
filtlist_short = sorted(list(dict.fromkeys(dict_filter_short)))
filtlist_long = sorted(list(dict.fromkeys(dict_filter_long)))
if len(dict_images_combined[d][f]) == 0:
dict_images_combined[d][f] = {'images': [file]}
else:
dict_images_combined[d][f]['images'].append(file)
print("Available Detectors for SW channel:", detlist_short)
print("Available Detectors for LW channel:", detlist_long)
print("Available SW Filters:", filtlist_short)
print("Available LW Filters:", filtlist_long)
Downloading images Available Detectors for SW channel: ['NRCB1'] Available Detectors for LW channel: ['NRCB5'] Available SW Filters: ['F070W', 'F115W', 'F200W'] Available LW Filters: ['F277W', 'F356W', 'F444W']
plt.figure(figsize=(14, 14))
for det in dets_short:
for i, filt in enumerate(filts_short):
image = fits.open(dict_images_combined[det][filt]['images'][0])
data_sb = image[1].data
ax = plt.subplot(1, len(filts_short), i + 1)
norm = simple_norm(data_sb, 'sqrt', percent=99.)
plt.xlabel("X [px]", fontdict=font2)
plt.ylabel("Y [px]", fontdict=font2)
plt.title(filt, fontdict=font2)
ax.imshow(data_sb, norm=norm, cmap='Greys')
plt.tight_layout()
dict_aper = {}
for det in dets_short:
dict_aper.setdefault(det, {})
for j, filt in enumerate(filts_short):
dict_aper[det].setdefault(filt, {})
dict_aper[det][filt]['stars for ap phot'] = None
dict_aper[det][filt]['stars for ap phot matched'] = None
dict_aper[det][filt]['aperture phot table'] = None
def find_bright_stars(det='NRCA1', filt='F070W', dist_sel=False):
bkgrms = MADStdBackgroundRMS()
mmm_bkg = MMMBackground()
image = fits.open(dict_images_combined[det][filt]['images'][i])
data_sb = image[1].data
imh = image[1].header
print("Selecting stars for aperture photometry on image {number} of filter {f}, detector {d}".format(number=i + 1, f=filt, d=det))
data = data_sb / imh['PHOTMJSR']
print("Conversion factor from {units} to DN/s for filter {f}:".format(units=imh['BUNIT'], f=filt), imh['PHOTMJSR'])
sigma_psf = dict_utils[filt]['psf fwhm']
print("FWHM for the filter {f}:".format(f=filt), sigma_psf, "px")
std = bkgrms(data)
bkg = mmm_bkg(data)
daofind = DAOStarFinder(threshold=th[j] * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
apcorr_stars = daofind(data)
dict_aper[det][filt]['stars for ap phot'] = apcorr_stars
if dist_sel:
print("")
print("Calculating closest neigbhour distance")
d = []
daofind_tot = DAOStarFinder(threshold=10 * std + bkg, fwhm=sigma_psf, roundhi=1.0, roundlo=-1.0,
sharplo=0.30, sharphi=1.40)
stars_tot = daofind_tot(data)
x_tot = stars_tot['xcentroid']
y_tot = stars_tot['ycentroid']
for xx, yy in zip(apcorr_stars['xcentroid'], apcorr_stars['ycentroid']):
sep = []
dist = np.sqrt((x_tot - xx)**2 + (y_tot - yy)**2)
sep = np.sort(dist)[1:2][0]
d.append(sep)
apcorr_stars['min distance'] = d
mask_dist = (apcorr_stars['min distance'] > min_sep[j])
apcorr_stars = apcorr_stars[mask_dist]
dict_aper[det][filt]['stars for ap phot'] = apcorr_stars
print("Minimum distance required:", min_sep[j], "px")
print("")
print("Number of bright isolated sources found in the image for {f}:".format(f=filt), len(apcorr_stars))
print("-----------------------------------------------------")
print("")
else:
print("")
print("Number of bright sources found in the image for {f}:".format(f=filt), len(apcorr_stars))
print("--------------------------------------------")
print("")
return
tic = time.perf_counter()
th = [700, 500] # threshold level for the two filters (length must match number of filters analyzed)
min_sep = [10, 10] # minimum separation acceptable for zp stars from closest neighbour
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):
find_bright_stars(det=det, filt=filt, dist_sel=False)
toc = time.perf_counter()
print("Elapsed Time for finding stars for Aperture Photometry:", toc - tic)
Selecting stars for aperture photometry on image 1 of filter F115W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F115W: 3.821892261505127 FWHM for the filter F115W: 1.298 px Number of bright sources found in the image for F115W: 1302 -------------------------------------------- Selecting stars for aperture photometry on image 1 of filter F200W, detector NRCB1 Conversion factor from MJy/sr to DN/s for filter F200W: 2.564082860946655 FWHM for the filter F200W: 2.141 px Number of bright sources found in the image for F200W: 1483 -------------------------------------------- Elapsed Time for finding stars for Aperture Photometry: 7.768850879999718
for det in dets_short:
for j, filt in enumerate(filts_short):
for i in np.arange(0, len(dict_images_combined[det][filt]['images']), 1):
image = ImageModel(dict_images_combined[det][filt]['images'][i])
ra, dec = image.meta.wcs(dict_aper[det][filt]['stars for ap phot']['xcentroid'],
dict_aper[det][filt]['stars for ap phot']['ycentroid'])
radec = SkyCoord(ra, dec, unit='deg')
dict_aper[det][filt]['stars for ap phot']['radec'] = radec
idx_ap, d2d_ap, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot']['radec'],
dict_aper[det][filt2]['stars for ap phot']['radec'])
sep_constraint_ap = d2d_ap < max_sep
matched_apcorr_f115w = Table()
matched_apcorr_f200w = Table()
matched_apcorr_f115w = dict_aper[det][filt1]['stars for ap phot'][sep_constraint_ap]
matched_apcorr_f200w = dict_aper[det][filt2]['stars for ap phot'][idx_ap[sep_constraint_ap]]
dict_aper[det][filt1]['stars for ap phot matched'] = matched_apcorr_f115w
dict_aper[det][filt2]['stars for ap phot matched'] = matched_apcorr_f200w
if os.path.isfile('./aperture_correction_table.txt'):
ap_tab = './aperture_correction_table.txt'
else:
print("Downloading the aperture correction table")
boxlink_apcorr_table = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/aperture_correction_table.txt'
boxfile_apcorr_table = './aperture_correction_table.txt'
urllib.request.urlretrieve(boxlink_apcorr_table, boxfile_apcorr_table)
ap_tab = './aperture_correction_table.txt'
aper_table = pd.read_csv(ap_tab, header=None, sep='\s+', index_col=0,
names=['filter', 'pupil', 'wave', 'r10', 'r20', 'r30', 'r40', 'r50', 'r60', 'r70', 'r80',
'r85', 'r90', 'sky_flux_px', 'apcorr10', 'apcorr20', 'apcorr30', 'apcorr40',
'apcorr50', 'apcorr60', 'apcorr70', 'apcorr80', 'apcorr85', 'apcorr90', 'sky_in',
'sky_out'], comment='#', skiprows=0, usecols=range(0, 26))
aper_table.head()
Downloading the aperture correction table
| pupil | wave | r10 | r20 | r30 | r40 | r50 | r60 | r70 | r80 | ... | apcorr30 | apcorr40 | apcorr50 | apcorr60 | apcorr70 | apcorr80 | apcorr85 | apcorr90 | sky_in | sky_out | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| filter | |||||||||||||||||||||
| F070W | CLEAR | 70 | 0.451 | 0.869 | 1.263 | 1.648 | 2.191 | 3.266 | 5.176 | 7.292 | ... | 3.3651 | 2.5305 | 2.0347 | 1.7210 | 1.5328 | 1.4174 | 1.4174 | 1.5880 | 7.292 | 9.017 |
| F090W | CLEAR | 90 | 0.408 | 0.638 | 0.992 | 1.503 | 1.925 | 2.549 | 4.162 | 7.480 | ... | 3.3519 | 2.5241 | 2.0253 | 1.6977 | 1.4908 | 1.4173 | 1.4173 | 1.6016 | 7.480 | 9.251 |
| F115W | CLEAR | 115 | 0.374 | 0.571 | 0.778 | 1.036 | 1.768 | 2.324 | 3.287 | 6.829 | ... | 3.3404 | 2.5070 | 2.0131 | 1.6825 | 1.4520 | 1.3310 | 1.3310 | 1.3964 | 6.829 | 9.723 |
| F140M | CLEAR | 140 | 0.414 | 0.617 | 0.801 | 1.031 | 1.367 | 2.434 | 3.118 | 6.031 | ... | 3.3397 | 2.5060 | 2.0067 | 1.6815 | 1.4465 | 1.3029 | 1.3029 | 1.3794 | 6.031 | 9.608 |
| F150W | CLEAR | 150 | 0.431 | 0.639 | 0.826 | 1.065 | 1.360 | 2.476 | 3.199 | 6.082 | ... | 3.3405 | 2.5067 | 2.0070 | 1.6828 | 1.4485 | 1.3068 | 1.3068 | 1.4188 | 6.082 | 9.496 |
5 rows × 25 columns
def aperture_phot(det=det, filt='F070W'):
radii = [aper_table.loc[filt]['r70']]
ees = '70'.split()
ee_radii = dict(zip(ees, radii))
positions = np.transpose((dict_aper[det][filt]['stars for ap phot matched']['xcentroid'],
dict_aper[det][filt]['stars for ap phot matched']['ycentroid']))
image = fits.open(dict_images_combined[det][filt]['images'][0])
data_sb = image[1].data
imh = image[1].header
data = data_sb / imh['PHOTMJSR']
# sky from the aperture correction table:
sky = {"sky_in": aper_table.loc[filt]['r80'], "sky_out": aper_table.loc[filt]['r85']}
tic = time.perf_counter()
table_aper = Table()
for ee, radius in ee_radii.items():
print("Performing aperture photometry for radius equivalent to EE = {0}% for filter {1}".format(ee, filt))
aperture = CircularAperture(positions, r=radius)
annulus_aperture = CircularAnnulus(positions, r_in=sky["sky_in"], r_out=sky["sky_out"])
annulus_mask = annulus_aperture.to_mask(method='center')
bkg_median = []
for mask in annulus_mask:
annulus_data = mask.multiply(data)
annulus_data_1d = annulus_data[mask.data > 0]
_, median_sigclip, _ = sigma_clipped_stats(annulus_data_1d)
bkg_median.append(median_sigclip)
bkg_median = np.array(bkg_median)
phot = aperture_photometry(data, aperture, method='exact')
phot['annulus_median'] = bkg_median
phot['aper_bkg'] = bkg_median * aperture.area
phot['aper_sum_bkgsub'] = phot['aperture_sum'] - phot['aper_bkg']
apcorr = [aper_table.loc[filt]['apcorr70']]
phot['aper_sum_corrected'] = phot['aper_sum_bkgsub'] * apcorr
phot['mag_corrected'] = -2.5 * np.log10(phot['aper_sum_corrected']) + dict_utils[filt]['VegaMAG zp modB']
table_aper.add_column(phot['aperture_sum'], name='aper_sum_' + ee)
table_aper.add_column(phot['annulus_median'], name='annulus_median_' + ee)
table_aper.add_column(phot['aper_bkg'], name='aper_bkg_ee_' + ee)
table_aper.add_column(phot['aper_sum_bkgsub'], name='aper_sum_bkgsub_' + ee)
table_aper.add_column(phot['aper_sum_corrected'], name='aper_sum_corrected_' + filt)
table_aper.add_column(phot['mag_corrected'], name='mag_corrected_' + filt)
dict_aper[det][filt]['aperture phot table'] = table_aper
toc = time.perf_counter()
print("Time Elapsed:", toc - tic)
return
aperture_phot(det=det, filt=filt1)
aperture_phot(det=det, filt=filt2)
Performing aperture photometry for radius equivalent to EE = 70% for filter F115W Time Elapsed: 0.5989792430000307 Performing aperture photometry for radius equivalent to EE = 70% for filter F200W Time Elapsed: 0.5241444849998516
plt.figure(figsize=(14, 8))
plt.clf()
ax1 = plt.subplot(2, 1, 1)
ax1.set_xlabel(filt1, fontdict=font2)
ax1.set_ylabel('Zeropoint', fontdict=font2)
idx_zp_1, d2d_zp_1, _ = match_coordinates_sky(dict_aper[det][filt1]['stars for ap phot matched']['radec'], radec_f115w_inst)
sep_constraint_zp_1 = d2d_zp_1 < max_sep
f115w_ap_matched = np.array(dict_aper[det][filt1]['aperture phot table']['mag_corrected_' + filt1][sep_constraint_zp_1])
f115w_psf_matched = np.array(df_f115w[filt1 + '_inst'][idx_zp_1[sep_constraint_zp_1]])
diff_f115w = f115w_ap_matched - f115w_psf_matched
_, zp_f115w, zp_sigma_f115w = sigma_clipped_stats(diff_f115w)
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f115w) - 0.5
ylim1 = np.mean(diff_f115w) + 0.5
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(f115w_psf_matched, diff_f115w, s=50, color='k')
ax1.plot([xlim0, xlim1], [zp_f115w, zp_f115w], color='r', lw=5, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + ' Zeropoint = %5.3f $\pm$ %5.3f' % (zp_f115w, zp_sigma_f115w), color='k', fontdict=font2)
ax2 = plt.subplot(2, 1, 2)
ax2.set_xlabel(filt2, fontdict=font2)
ax2.set_ylabel('Zeropoint', fontdict=font2)
idx_zp_2, d2d_zp_2, _ = match_coordinates_sky(dict_aper[det][filt2]['stars for ap phot matched']['radec'], radec_f200w_inst)
sep_constraint_zp_2 = d2d_zp_2 < max_sep
f200w_ap_matched = np.array(dict_aper[det][filt2]['aperture phot table']['mag_corrected_' + filt2][sep_constraint_zp_2])
f200w_psf_matched = np.array(df_f200w[filt2 + '_inst'][idx_zp_2[sep_constraint_zp_2]])
diff_f200w = f200w_ap_matched - f200w_psf_matched
_, zp_f200w, zp_sigma_f200w = sigma_clipped_stats(diff_f200w)
xlim0 = -9
xlim1 = -5
ylim0 = np.mean(diff_f200w) - 0.5
ylim1 = np.mean(diff_f200w) + 0.5
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(f200w_psf_matched, diff_f200w, s=50, color='k')
ax2.plot([xlim0, xlim1], [zp_f200w, zp_f200w], color='r', lw=5, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + ' Zeropoint = %5.3f $\pm$ %5.3f' % (zp_f200w, zp_sigma_f200w), color='k', fontdict=font2)
plt.tight_layout()
if os.path.isfile('./pointsource.cat'):
input_cat = './pointsource.cat'
else:
print("Downloading input pointsource catalog")
boxlink_input_cat = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/pointsource.cat'
boxfile_input_cat = './pointsource.cat'
urllib.request.urlretrieve(boxlink_input_cat, boxfile_input_cat)
input_cat = './pointsource.cat'
cat = pd.read_csv(input_cat, header=None, sep='\s+', names=['ra_in', 'dec_in', 'f070w_in', 'f115w_in',
'f200w_in', 'f277w_in', 'f356w_in', 'f444w_in'],
comment='#', skiprows=7, usecols=range(0, 8))
cat.head()
Downloading input pointsource catalog
| ra_in | dec_in | f070w_in | f115w_in | f200w_in | f277w_in | f356w_in | f444w_in | |
|---|---|---|---|---|---|---|---|---|
| 0 | 80.386396 | -69.468909 | 21.34469 | 20.75333 | 20.25038 | 20.23116 | 20.20482 | 20.23520 |
| 1 | 80.385588 | -69.469201 | 20.12613 | 19.33709 | 18.64676 | 18.63521 | 18.58796 | 18.64291 |
| 2 | 80.380365 | -69.470930 | 21.52160 | 20.98518 | 20.53500 | 20.51410 | 20.49231 | 20.51610 |
| 3 | 80.388130 | -69.468453 | 20.82162 | 20.06542 | 19.40552 | 19.39262 | 19.34927 | 19.40018 |
| 4 | 80.388936 | -69.468196 | 21.47197 | 20.92519 | 20.46507 | 20.44447 | 20.42186 | 20.44687 |
lim_ra_min = np.min(radec_f115w.ra)
lim_ra_max = np.max(radec_f115w.ra)
lim_dec_min = np.min(radec_f115w.dec)
lim_dec_max = np.max(radec_f115w.dec)
cat_sel = cat[(cat['ra_in'] > lim_ra_min) & (cat['ra_in'] < lim_ra_max) & (cat['dec_in'] > lim_dec_min)
& (cat['dec_in'] < lim_dec_max)]
plt.figure(figsize=(12, 14))
plt.clf()
ax1 = plt.subplot(1, 2, 1)
mag1_in = np.array(cat_sel['f115w_in'])
mag2_in = np.array(cat_sel['f200w_in'])
diff_in = mag1_in - mag2_in
xlim0 = -0.25
xlim1 = 1.75
ylim0 = 25
ylim1 = 15
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(mag1_in - mag2_in, mag1_in, s=1, color='k')
ax1.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)
ax1.set_ylabel(filt1, fontdict=font2)
ax1.text(xlim0 + 0.15, 15.5, "Input", fontdict=font2)
ax2 = plt.subplot(1, 2, 2)
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
f115w = f115w_inst + zp_f115w
f200w = f200w_inst + zp_f200w
maglim = np.arange(18, 25, 1)
mags = []
errs_mag = []
errs_col = []
for i in np.arange(0, len(maglim) - 1, 1):
mag = (maglim[i] + maglim[i + 1]) / 2
err_mag1 = ef115w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]
err_mag2 = ef200w_inst[(f115w > maglim[i]) & (f115w < maglim[i + 1])]
err_mag = np.mean(err_mag1[i])
err_temp = np.sqrt(err_mag1**2 + err_mag2**2)
err_col = np.mean(err_temp[i])
errs_mag.append(err_mag)
errs_col.append(err_col)
mags.append(mag)
col = [0] * (len(maglim) - 1)
ax2.errorbar(col, mags, yerr=errs_mag, xerr=errs_col, fmt='o', color='k')
ax2.scatter(f115w - f200w, f115w, s=1, color='k')
ax2.text(xlim0 + 0.15, 15.5, "Output", fontdict=font2)
ax2.set_xlabel(filt1 + ' - ' + filt2, fontdict=font2)
ax2.set_ylabel(filt1, fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(14, 8))
plt.clf()
ax1 = plt.subplot(2, 1, 1)
ax1.set_xlabel(filt1, fontdict=font2)
ax1.set_ylabel('$\Delta$ Mag', fontdict=font2)
radec_input = SkyCoord(cat_sel['ra_in'], cat_sel['dec_in'], unit='deg')
idx_f115w_cfr, d2d_f115w_cfr, _ = match_coordinates_sky(radec_input, radec_f115w)
sep_f115w_cfr = d2d_f115w_cfr < max_sep
f115w_inp_cfr = np.array(cat_sel['f115w_in'][sep_f115w_cfr])
f115w_psf_cfr = np.array(f115w[idx_f115w_cfr[sep_f115w_cfr]])
diff_f115w_cfr = f115w_inp_cfr - f115w_psf_cfr
_, med_diff_f115w_cfr, sig_diff_f115w_cfr = sigma_clipped_stats(diff_f115w_cfr)
xlim0 = 16
xlim1 = 24.5
ylim0 = np.mean(diff_f115w_cfr) - 0.5
ylim1 = np.mean(diff_f115w_cfr) + 0.5
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.scatter(f115w_psf_cfr, diff_f115w_cfr, s=5, color='k')
ax1.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 0.15, filt1 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f'
% (med_diff_f115w_cfr, sig_diff_f115w_cfr), color='k', fontdict=font2)
ax2 = plt.subplot(2, 1, 2)
ax2.set_xlabel(filt2, fontdict=font2)
ax2.set_ylabel('$\Delta$ Mag', fontdict=font2)
idx_f200w_cfr, d2d_f200w_cfr, _ = match_coordinates_sky(radec_input, radec_f200w)
sep_f200w_cfr = d2d_f200w_cfr < max_sep
f200w_inp_cfr = np.array(cat_sel['f200w_in'][sep_f200w_cfr])
f200w_psf_cfr = np.array(f200w[idx_f200w_cfr[sep_f200w_cfr]])
diff_f200w_cfr = f200w_inp_cfr - f200w_psf_cfr
_, med_diff_f200w_cfr, sig_diff_f200w_cfr = sigma_clipped_stats(diff_f200w_cfr)
xlim0 = 16
xlim1 = 24
ylim0 = np.mean(diff_f200w_cfr) - 0.5
ylim1 = np.mean(diff_f200w_cfr) + 0.5
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.scatter(f200w_psf_cfr, diff_f200w_cfr, s=5, color='k')
ax2.plot([xlim0, xlim1], [0, 0], color='r', lw=5, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 0.15, filt2 + ' $\Delta$ Mag = %5.3f $\pm$ %5.3f'
% (med_diff_f200w_cfr, sig_diff_f200w_cfr), color='k', fontdict=font2)
plt.tight_layout()
plt.figure(figsize=(12, 6))
ax1 = plt.subplot(1, 2, 1)
xlim0 = -10
xlim1 = 10
ylim0 = -10
ylim1 = 10
ax1.set_xlim(xlim0, xlim1)
ax1.set_ylim(ylim0, ylim1)
ax1.xaxis.set_major_locator(ticker.AutoLocator())
ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.yaxis.set_major_locator(ticker.AutoLocator())
ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax1.set_xlabel('$\Delta$ RA (mas)', fontdict=font2)
ax1.set_ylabel('$\Delta$ Dec (mas)', fontdict=font2)
ax1.set_title(filt1, fontdict=font2)
ra_f115w_inp_cfr = np.array(cat_sel['ra_in'][sep_f115w_cfr])
ra_f115w_psf_cfr = np.array(radec_f115w.ra[idx_f115w_cfr[sep_f115w_cfr]])
dec_f115w_inp_cfr = np.array(cat_sel['dec_in'][sep_f115w_cfr])
dec_f115w_psf_cfr = np.array(radec_f115w.dec[idx_f115w_cfr[sep_f115w_cfr]])
dec_rad_f115w = np.radians(dec_f115w_psf_cfr)
diffra_f115w_cfr = ((((ra_f115w_inp_cfr - ra_f115w_psf_cfr) * np.cos(dec_rad_f115w)) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffra_f115w_cfr, sig_diffra_f115w_cfr = sigma_clipped_stats(diffra_f115w_cfr)
diffdec_f115w_cfr = (((dec_f115w_inp_cfr - dec_f115w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffdec_f115w_cfr, sig_diffdec_f115w_cfr = sigma_clipped_stats(diffdec_f115w_cfr)
ax1.scatter(diffra_f115w_cfr, diffdec_f115w_cfr, s=1, color='k')
ax1.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax1.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax1.text(xlim0 + 0.05, ylim1 - 1.50, ' $\Delta$ RA (mas) = %5.3f $\pm$ %5.3f'
% (med_diffra_f115w_cfr, sig_diffra_f115w_cfr), color='k', fontdict=font2)
ax1.text(xlim0 + 0.05, ylim1 - 3.0, ' $\Delta$ Dec (mas) = %5.3f $\pm$ %5.3f'
% (med_diffdec_f115w_cfr, sig_diffdec_f115w_cfr), color='k', fontdict=font2)
ax2 = plt.subplot(1, 2, 2)
xlim0 = -10
xlim1 = 10
ylim0 = -10
ylim1 = 10
ax2.set_xlim(xlim0, xlim1)
ax2.set_ylim(ylim0, ylim1)
ax2.set_title(filt2, fontdict=font2)
ax2.xaxis.set_major_locator(ticker.AutoLocator())
ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.yaxis.set_major_locator(ticker.AutoLocator())
ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())
ax2.set_xlabel('$\Delta$ RA (mas)', fontdict=font2)
ax2.set_ylabel('$\Delta$ Dec (mas)', fontdict=font2)
ra_f200w_inp_cfr = np.array(cat_sel['ra_in'][sep_f200w_cfr])
ra_f200w_psf_cfr = np.array(radec_f200w.ra[idx_f200w_cfr[sep_f200w_cfr]])
dec_f200w_inp_cfr = np.array(cat_sel['dec_in'][sep_f200w_cfr])
dec_f200w_psf_cfr = np.array(radec_f200w.dec[idx_f200w_cfr[sep_f200w_cfr]])
dec_rad_f200w = np.radians(dec_f200w_psf_cfr)
diffra_f200w_cfr = ((((ra_f200w_inp_cfr - ra_f200w_psf_cfr) * np.cos(dec_rad_f200w)) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffra_f200w_cfr, sig_diffra_f200w_cfr = sigma_clipped_stats(diffra_f200w_cfr)
diffdec_f200w_cfr = (((dec_f200w_inp_cfr - dec_f200w_psf_cfr) * u.deg).to(u.mas) / (1 * u.mas))
_, med_diffdec_f200w_cfr, sig_diffdec_f200w_cfr = sigma_clipped_stats(diffdec_f200w_cfr)
ax2.scatter(diffra_f200w_cfr, diffdec_f200w_cfr, s=1, color='k')
ax2.plot([0, 0], [ylim0, ylim1], color='k', lw=2, ls='--')
ax2.plot([xlim0, xlim1], [0, 0], color='k', lw=2, ls='--')
ax2.text(xlim0 + 0.05, ylim1 - 1.50, ' $\Delta$ RA (mas) = %5.3f $\pm$ %5.3f'
% (med_diffra_f200w_cfr, sig_diffra_f200w_cfr), color='k', fontdict=font2)
ax2.text(xlim0 + 0.05, ylim1 - 3.0, ' $\Delta$ Dec (mas) = %5.3f $\pm$ %5.3f'
% (med_diffdec_f200w_cfr, sig_diffdec_f200w_cfr), color='k', fontdict=font2)
plt.tight_layout()